CN104237871B - Delay inequality estimation method based on phase compensation - Google Patents

Abstract

The invention relates to a delay inequality estimation method based on phase compensation. The method is used for estimating the delay inequality of different receiving array elements for receiving target radiation signals in an underwater sound ternary array passive positioning system. The method includes the steps that phase compensation is performed on the cross-power spectrum signals between the receiving array elements; energy obtaining and peak search are performed on the cross-power spectrum signals subjected to phase compensation, and a delay inequality period value is obtained; the delay inequality correction value of the cross-power spectrum signals subjected to phase compensation is obtained on the phase compensation peak point through a least square method; the delay inequality period value and the delay inequality correction value are added together, and a final delay inequality estimation value is obtained.

Description

Delay inequality estimation method based on phase compensation

Technical Field

The invention relates to the field of sonar signal processing, in particular to a time delay difference estimation method based on phase compensation.

Background

The use of radiated noise from underwater targets to detect and locate the targets is an important means of submarine covert discovery and attack. At present, the passive positioning method mainly comprises 3 types of ternary arrays, target motion analysis and matching field processing, wherein the ternary passive positioning method is still a passive positioning method generally adopted on submarine sonar equipment, and the key is two time delay difference estimations between 3 sub-array received signals, because the time delay difference estimation accuracy directly determines the positioning accuracy of a target.

Because most of the radiated noise of the underwater target is broadband signals, the cross-correlation function of the signals generally has a main lobe with a sharp correlation peak, and the receiving end can obtain the time delay difference between the target radiated signal and each array element by performing cross-correlation on the received signals. In order to break through the limitation of sampling rate, researchers have proposed a variety of methods for refining correlation peaks, which are mainly classified into two types: one is that the local thinning method is directly carried out on the peak value of the correlation function, such as cubic spline interpolation method, quadratic spline interpolation method, cosine function interpolation method, etc., and the method has small operand but low precision; the other type is that the correlation peak is refined by preprocessing, such as an MCZT method, a ZOOM-FFT method, a BCZT method and the like, the method does not interpolate the peak locally, but realizes the refinement by frequency domain operation, the operation amount is large, but the precision is high.

In order to improve the estimation accuracy of the time delay difference between two paths of signals, J.S. Berndatt and A.G. Pielshol propose to estimate the time delay difference between the signals by calculating the cross power spectrum of the two signals, but the phase ambiguity problem occurs in practical application, the tolerance is poor, and the estimation accuracy of the time delay difference is not high.

Disclosure of Invention

The invention aims to solve the problems of phase ambiguity, poor tolerance, low delay inequality estimation precision and the like in practical application of the delay inequality estimation method in the prior art, thereby providing the delay inequality estimation method with high tolerance and high delay inequality estimation precision.

In order to achieve the above object, the present invention provides a delay inequality estimation method based on phase compensation, which is used for estimating delay inequality of target radiation signals received by different receiving array elements in an underwater acoustic ternary array passive positioning system, and the method includes:

performing phase compensation on the cross-power spectrum signals among the receiving array elements on a frequency domain;

then, energy calculation and peak value search are carried out on the cross-power spectrum signal after phase compensation, and a time delay difference period value is obtained;

then, a time delay difference correction value is obtained for the cross-spectrum signal after the phase compensation according to a least square method at the peak point of the phase compensation;

and finally, adding the delay difference period value and the delay difference correction value to obtain a final delay difference estimation value.

In the above technical solution, the method further comprises:

step 1), receiving a target radiation signal by each array element in a receiving array, and then preprocessing the target radiation signal to obtain a frequency domain signal of the target radiation signal received by each array element;

step 2), selecting any two receiving array elements, and multiplying frequency domain signals of target radiation signals received by the two receiving array elements after Fourier transformation to obtain frequency domain cross spectrum data of the two receiving array elements;

step 3), performing phase compensation on the frequency domain cross-spectrum data obtained in the step 2);

step 4), energy is obtained for the result obtained in the step 3) after phase compensation, then energy peak value search is carried out according to the result obtained by energy peak value search, and the time value corresponding to the obtained energy peak value point is a delay difference period value;

step 5), solving a phase difference of the frequency domain cross spectrum data after the phase compensation at the energy peak point obtained after the phase compensation, and then solving a time delay difference correction value of the phase difference according to a least square method;

and 6) adding the period value of the time delay difference obtained in the step 4) and the time delay difference correction value obtained in the step 5) to obtain the final time delay difference estimation value difference of the two receiving array elements selected in the step 2).

In the above technical solution, in the step 1), the pretreatment includes:

step 1-1), filtering and amplifying target radiation signals received by each array element;

step 1-2), carrying out discrete sampling on the filtered signal at a certain sampling rate;

and 1-3) carrying out Fourier transform on the discrete sampling result to obtain a corresponding frequency domain signal.

In the above technical solution, in the step 3), a calculation formula for performing phase compensation on the frequency domain cross spectrum data obtained in the step 2) is as follows:

${G^{\′}}_{x_1{x}_2}(\mathrm{\Δt},f)=G_{x_1{x}_2}\left(f\right)e^{j2\mathrm{\π\Δtf}},f_l\≤f\≤f_h,-\frac{d}{c}\≤\mathrm{\Δt}=\frac{1}{{\mathrm{kf}}_s}\≤\frac{d}{c}$

sampling rate improvement factor, f, in phase compensation_lLower limit of filter, f_hIs the upper limit of the filter, f_sIs the sampling rate.

The method enhances the compatibility of the cross-spectrum method in practical application and improves the estimation precision of the time delay difference.

The invention has the advantages that:

the method can obtain the time delay difference with high precision, has strong practicability, and can effectively improve the positioning precision of the underwater acoustic ternary array passive positioning system.

Drawings

FIG. 1 is a schematic diagram of a ternary array involved in the present invention;

FIG. 2 is a flow chart of the method of the present invention;

FIG. 3 is a diagram of the standard deviation of delay variation estimation errors for different methods;

fig. 4(a) is the estimation result of delay inequality between array elements 1 and 2 when the target and the array are static;

fig. 4(b) is the estimation result of delay inequality between array elements 2 and 3 when the target and the base array are static;

fig. 5(a) is the estimation result of the delay inequality between the array elements 1 and 2 when the target and the matrix move relatively;

fig. 5(b) shows the estimation result of the delay inequality between the array elements 2 and 3 when the target and the matrix move relatively.

Detailed Description

The invention will now be further described with reference to the accompanying drawings.

Before describing the method of the present invention in detail, a receiving array to which the method of the present invention is applied will be described first.

Fig. 1 is a schematic diagram of a ternary array, which includes three array elements, denoted by numerals 1, 2, and 3, respectively, and a target radiates a signal from the θ direction, and the radiated signal reaches each array element after propagating through an underwater acoustic channel. For convenience of explanation, the method of the present invention will be described below by taking 1 and 2 array elements as examples, and the method of the present invention is also applicable to other array elements.

Based on a receiving array such as that shown in fig. 1, the invention firstly directly performs phase compensation on cross-power spectrum signals among array elements in a frequency domain; then, energy calculation and peak value search are carried out on the compensated cross-power spectrum signal to obtain a time delay difference period value; then, a time delay difference correction value is obtained for the compensated cross-spectrum signal at the peak point of the phase compensation according to a least square method; and finally, adding the delay difference correction values of the delay difference period value domain to obtain a final delay difference estimation value.

Referring to fig. 2, the method of the present invention specifically includes:

step 1), receiving target radiation signals by each array element in the receiving array, and then preprocessing the target radiation signals to obtain frequency domain signals of the target radiation signals received by each array element.

The preprocessing of the target radiation signal includes:

step 1-1), filtering and amplifying target radiation signals received by each array element, and marking the filtered signals as x_m(t), wherein m represents a matrix element number;

step 1-2), at a sampling rate f_sFor the filtered signal x_m(t) discrete sampling to obtain x_m(nT_s) (ii) a Wherein T is_s＝1f_sRepresenting a sampling time interval;

step 1-3), discrete sampling result x_m(nT_s) Fourier transform is carried out to obtain corresponding frequency domain signal X_m(f),f_l≤f≤f_h(ii) a Wherein f is_lLower limit of filter, f_hIs the filter upper bound;

step 2), Fourier of target radiation signals received by the receiving array element 1 and the receiving array element 2Transformed frequency domain signal X₁(f),f_l≤f≤f_hAndf_l≤f≤f_hmultiplying to obtain frequency domain cross spectrum data of the two receiving array elementsf_l≤f≤f_hWherein^*() Is complex conjugation.

Step 3) of the frequency domain cross spectrum data obtained in the step 2)f_l≤f≤f_hAnd carrying out phase compensation.

The calculation formula for performing phase compensation on the frequency domain cross-spectrum data is as follows:

${G^{\′}}_{x_1{x}_2}(\mathrm{\Δt},f)=G_{x_1{x}_2}\left(f\right)e^{j2\mathrm{\π\Δtf}},f_l\≤f\≤f_h,-\frac{d}{c}\≤\mathrm{\Δt}=\frac{1}{{\mathrm{kf}}_s}\≤\frac{d}{c}---\left(1\right)$

wherein, ${G^{\′}}_{x_1{x}_2}(\mathrm{\Δt},f),f_l\≤f\≤f_h,-\frac{d}{c}\≤\mathrm{\Δt}=\frac{1}{{\mathrm{kf}}_s}\≤\frac{d}{c}$ and (4) representing a result of phase compensation of the frequency domain cross-spectrum data, wherein d is the array element interval, c is the effective sound velocity, and k is the sampling rate improvement multiple during phase compensation.

And 4) performing energy calculation on the result obtained in the step 3) after the phase compensation, and then performing energy peak value search according to the result of the energy calculation, wherein the time value corresponding to the obtained energy peak value point is a delay difference period value.

In this step, the formula for energy calculation of the result after phase compensation is shown in formula (2):

${{\mathrm{Pow}}^{\′}}_{12}\left(\mathrm{\Δt}\right)=\underset{f=f_l}{\overset{f_h}{\mathrm{\Σ}}}\left|{X^{\′}}_{12}(\mathrm{\Δt},f)\right|,-\frac{d}{c}\≤\mathrm{\Δt}=\frac{1}{{\mathrm{kf}}_s}\≤\frac{d}{c}---\left(2\right)$

the calculation formula of the energy peak search is shown in formula (3):

$\left\{\begin{array}{c}{{\mathrm{Pow}}^{\′}}_{12\max }\left(\mathrm{\Δtindex}\right)=\max \left({{\mathrm{Pow}}^{\′}}_{12}\left(\mathrm{\Δt}\right)\right),-\frac{d}{c}\≤\mathrm{\Δt}=\frac{1}{{\mathrm{kf}}_s}\≤\frac{d}{c}\\ {\mathrm{\τ}}_{12N}=\mathrm{\Δtindex},\end{array}\right.---\left(3\right)$

wherein tau is_12NIndicating the delay difference period value.

And 5) solving a phase difference of the frequency domain cross spectrum data after the phase compensation at the energy peak point obtained after the phase compensation, and then solving a time delay difference correction value of the phase difference according to a least square method.

In this step, the formula for obtaining the phase difference is shown in the following formula (4):

thereinf_l≤f≤f_hIndicating the phase difference.

Wherein tau is_12MThe delay difference correction value is represented.

Step 6), the period value tau of the time delay difference obtained in the step 4) is used for calculating_12NThe corrected value tau of the time delay difference obtained in step 5)_12MAdding to obtain the final time delay difference estimation value difference

${\widehat{\mathrm{\τ}}}_{12}={\mathrm{\τ}}_{12N}+{\mathrm{\τ}}_{12M}---\left(6\right)$

The above is a description of the basic steps of the method of the present invention, and the necessity of the delay inequality estimation method of the present invention is explained below.

Assuming that a target radiation signal is a band-limited signal s (T), the length of the signal is T, and the upper and lower frequency limits of the signal are f_lAnd f_h. LetterThe number arriving at the receiving end is:

$\left\{\begin{array}{c}{x}_1\left(t\right)=k_{1}s(t-{\mathrm{\τ}}_1)+n_1\left(t\right),\\ x_2\left(t\right)=k_{2}s(t-{\mathrm{\τ}}_2)+n_2\left(t\right).\end{array}\right.---\left(7\right)$

wherein x is₁(t) and x₂(t) receiving signals, tau, for array elements 1, 2, respectively₁And τ₂Propagation delay, k, of the target radiation signal to array elements 1, 2, respectively₁And k₂Amplitude, n, of two signals respectively₁(t) and n₂(t) are white Gaussian noises, respectively.

Firstly, Fourier transform is carried out on two paths of signals to obtain:

$\left\{\begin{array}{c}{X}_1\left(f\right)={\∫}_0^T{x}_1\left(t\right)e^{j2\mathrm{\πft}}\mathrm{dt}=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{x}_1\left(n\right)e^{j2\mathrm{\πfn}},\\ X_2\left(f\right)={\∫}_0^T{x}_2\left(t\right)e^{j2\mathrm{\πft}}\mathrm{dt}=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{x}_2\left(n\right)e^{j2\mathrm{\πfn}}.\end{array}\right.1\≤f\≤f_s,N=T\·f_s---\left(8\right)$

further, the cross power spectrum between the two paths of signals can be obtained as follows:

$G_{x_1{x}_2}\left(f\right)=X_1\left(f\right)\·X_2^{*}\left(f\right),f_1\≤f\≤f_h---\left(9\right)$

wherein (C)^*Is a conjugate transpose. The phase of the cross-power spectrum signal obtained according to equation (9) is:

from formula (10): phase positionAlready containing the delay difference tau of the received signals of the array elements 1, 2₁₂To, forUsing least squaresFitting is performed to determine τ₁₂The optimal estimate of (c) is:

if it is directly adoptedAs a final delay difference estimation, when the distance between the array elements 1 and 2 is more than half a wavelength, this will result inA phase ambiguity condition occurs, in which case it cannot be directly caused byThe periodic value (main value) of the phase difference between the signals received by the array elements 1 and 2 can not be measured directlyMeasuring the real time delay difference of the signals received by the array elements 1 and 2; aiming at the problem, the cross power spectrum of the signals received by the phase compensation array elements 1 and 2 can be adopted to solve the problem, then the energy calculation and peak search are carried out on the compensated cross power spectrum, and the delay difference period value (main value) of the signals received by the array elements 1 and 2 can be obtained, and the method comprises the following steps:

the phase compensation formula can be expressed by equation (12),

$G_{x_1{x}_2}^{\′}(t,f)=G_{x_1{x}_2}\left(f\right)e^{j2\mathrm{\πft}},f_1\≤f\≤f_h,t\∈[-\frac{d}{c}\≤\frac{1}{{\mathrm{kf}}_s}\≤\frac{d}{c}]---\left(12\right)$

wherein T ∈ [ -T': 1/(k.f)_s):T']For the phase compensation time, T' for the compensation cut-off time, f_sIs the sampling rate, 1/(k.f)_s) To compensate for the time interval, k in the equation is similar to a frequency interpolation method in order to improve the measurement of the delay difference period with high precision.

Then pair at each time pointf_l≤f≤f_hEnergy is obtained to obtain a set of energy values Pow (T), T ∈ [ -T ', T']Searching pow (t) to measure the period value tau of delay difference between array elements 1 and 2_12N. The compensation is limited by the sampling rate, which results in the period value tau of the delay difference_12NAnd the deviation exists between the real value of the delay difference, and the deviation is defined as the corrected value of the delay difference,

τ_12M＝τ₁₂-τ_12N, (13)

at this time, τ_12NThe value of (a) already ensures the calculation of Δ τ_12MWithout phase ambiguity, i.e.

Then, toAccording to the formula (11), the time delay difference correction value tau can be obtained by fitting by adopting a least square method_12MAnd (6) optimal estimation.

Finally, the time delay difference period value and the time delay difference correction value are added to obtain the time delay difference estimation value

τ₁₂＝τ_12N+τ_12M, (15)

As can be seen from the above derived equations, it is critical to realize the method to ensure that equation (14) holds, let τ be₁₂And τ_12NThe maximum allowable error is:

τ_12E＝max(|τ_12M|)＝max(|τ₁₂-τ_12N|), (16)

according to equation (14), the condition sufficient for no phase ambiguity is:

2πf_hτ_12E＜π, (17)

further simplification can result:

f_hτ_12E＜12, (18)

in practical application, the sampling rate generally satisfies Nyquest sampling theorem, and then the sampling rate f_s＞2f_hWhen this occurs, equation (18) can be rewritten as:

f_sτ_12E＜1, (19)

the physical significance of the formula (15) is that it is obtained if phase compensation is carried outThe detection error of the energy peak point is less than one sampling point, and the time delay difference correction value tau is calculated_12MUse ofThere is no phase ambiguity problem.

Sampling rate f in real systems_s＝Nf_hAnd in phase compensation equation (12), the sampling rate used for compensating the phase is generally relatively large (i.e., k ≧ 10), equation (19) can be rewritten as:

f_sτ_12E＜Nk2, (20)

equation (20) shows that as long as the error of the delay difference period value is less than the time length corresponding to Nk2 sampling points, the method can be used without phase ambiguity problem.

The implementation of the method of the present invention is described below with reference to an example.

Firstly, the received signal is filtered and amplified and is sampled at a sampling rate f_sSampling the filtered and amplified signal to obtain x₁(nT_s) And x₂(nT_s) Then, Fourier transform processing is performed on the sampled signal, and for the sake of calculation, at x₁(nT_s) And x₂(nT_s) And (3) complementing 0 at the tail part to enable the length of the two sequences to be N, wherein N should satisfy the minimum integer of the following conditions: 1) n is more than 2 (d.fs/c) +1, wherein d is the array spacing and c is the sound velocity; 2) log (log)₂N is an integer.

For x₁(nT_s) And x₂(nT_s) Performing fast fourier transform to:

$\left\{\begin{array}{c}{X}_1\left(m{w}_s\right)=\mathrm{FFT}\left(x_1\left(n{T}_s\right)\right),\\ X_2\left(m{w}_s\right)=\mathrm{FFT}\left(x_2\left(n{T}_s\right)\right),\end{array}\right.(0\&le;m<N)---\left(21\right)$

wherein, w_s＝2πNT_sIs the corresponding frequency sampling interval.

Is provided withThen by M₁To M₂The frequencies corresponding to the frequency points are all f_lTo f_hInternal;androunding up and down, respectively.

Secondly, the cross power spectrum of the two paths of signals is obtained according to the formula (9)M₁≤m≤M₂The cross-power spectrum is then phase compensated as per equation (22).

$G_{x_1{x}_2}^{\&prime;}(t,m{w}_s)=G_{x_1{x}_2}\left(m{w}_s\right)e^{j2\mathrm{\&pi;mt}},(M_1\&le;m\&le;M_2,-\frac{d}{c}\&le;t=\frac{1}{k\&CenterDot;\mathrm{fs}}\&le;\frac{d}{c})---\left(22\right)$

To pairEnergy is obtained at each time point, and then the array element 1 and 2 time delay difference period value tau is obtained by detecting the peak value_12N. Then processing the phase signal at the peak point compensation position according to the formula (10) and the formula (11) to obtain a time delay difference correction value tau_12M。

Finally, the delay difference period value tau is measured_12NSum delay difference correction value tau_12MAdding to obtain the final time delay difference estimated value tau₁₂。

In order to illustrate that the method of the invention is closer to the Cramer-Rao lower limit than the conventional method, and the time delay difference estimation precision is high, the Cramer-Rao lower limit of the time delay difference estimation standard difference calculated by the model is given below.

G.H. Konapus has pointed out that the Cramer-Rao lower limit of the delay inequality estimation standard deviation calculated by the model meets the following requirements:

$\mathrm{\&sigma;}\left(\mathrm{\&Delta;}\widetilde{\mathrm{\&tau;}}\right)=\{2T\underset{2\mathrm{\&pi;}{f}_l}{\overset{2\mathrm{\&pi;}{f}_h}{\&Integral;}}{w}^2\frac{|r\left(f\right){|}^2}{1-|r\left(f\right){|}^2}\mathrm{df}{\}}^{-1/2},---\left(23\right)$

wherein T is the sampling time, r (f) is x₁(t) and x₂(t) a coherence spectrum, defined as follows:

$|r\left(f\right){|}^2=\frac{K_{\mathrm{ss}}^2\left(f\right)}{[K_{\mathrm{ss}}\left(f\right)+K_{\mathrm{nn}}\left(f\right){]}^2},---\left(24\right)$

wherein, K_ss(f) And K_nn(f) Power spectra of the signal and noise, respectively; in the simple case, if K_ss(f) And K_nn(f) Are all flat in shape, i.e.

K_ss(f)＝S₀,K_nn(f)＝N₀,f∈(f_l,f_h) (25)

At this time, it is possible to obtain:

$\mathrm{\&sigma;}\left(\mathrm{\&Delta;}\widetilde{\mathrm{\&tau;}}\right)\&ap;[{\frac{3(N_0^2+2N_0{S}_0)}{8{\mathrm{\&pi;}}^{2}T{S}_0^2(f_h^3-f_l^3)}]}^{1/2},---\left(26\right)$

and the formula (26) is the Cramer-Rao lower limit of the estimation standard deviation of the passive sonar time delay difference.

Next, MATLAB is adopted to simulate the time delay difference estimation precision of different methods, broadband white Gaussian noise is selected as a target radiation signal in the simulation, and the upper and lower limits of the noise frequency are f_l500Hz and f_h1kHz, array element spacing d is 50m, target orientation theta is 60 degrees, target distance R is 2000m, effective sound velocity c is 1500m/s, and system early sampling rate f is_s10kHz, sample length T0.25 s, and after 2 down-sampling, the sampling rate is changed to f_s＝5kHz。

Fig. 3 shows the standard deviation of the delay difference estimation error in different methods, and it can be known from fig. 3 that:

1) the time delay difference estimation precision of the correlation method without interpolation is the lowest because of the limitation of the sampling rate;

2) compared with a correlation method without interpolation, the method for carrying out quadratic spline interpolation on the local part of the correlated peak envelope breaks through the limitation of sampling rate, the estimation precision of the time delay difference is greatly improved, but because much information is lost in the process of obtaining the peak envelope in the early stage, the method cannot achieve the high estimation precision of the time delay difference;

3) compared with a quadratic spline interpolation method, the ZOOM-FFT method related peak refining method is to refine related data in a frequency domain, so that the limitation that only data near a peak envelope is refined is avoided;

4) compared with other methods, the standard deviation of the time delay difference estimation error of the method is minimum and is closer to the Cramer-Rao lower limit, and particularly, the time delay difference estimation precision of the method is not obviously reduced under the condition of low signal-to-noise ratio; the main reason is that the method combines the advantages of the cross-spectrum method and the related peak refinement, and improves the combination process of the cross-spectrum method and the related peak refinement to a certain extent, and the method adopts related data to directly calculate the delay difference period value and the delay difference correction value in the frequency domain, and does not need to convert the related peak refinement method to the time domain to calculate the delay difference period value and the delay difference correction value, thereby reducing the operation amount, ensuring the precision and simultaneously ensuring that useful information is not lost. The simulation result shows that under the condition of selecting the parameters (the sampling rate is 10)³Hz), the time delay difference estimation precision of the method can reach 10^-6s。

In order to verify how the application effect of the method is in practical conditions, MATLAB is firstly adopted for numerical simulation analysis, broadband white Gaussian noise is selected as a target radiation signal in the following simulation under the simulation conditions, and the upper and lower limits of the noise frequency are f_l500Hz and f_h1kHz, array element spacing d is 50m, target orientation theta is 60 degrees, target distance R is 2000m, effective sound velocity c is 1500m/s, and system early-stage sampling rate is f_s10kHz, T0.25 s, snr 20dB, and f after 2 times down sampling_s5 kHz. According to the simulation analysis conditions, the following conditions are known: the time delay difference between the array elements 1 and 2 is a true value and is tau₁₂1631.0us, the time delay difference between the array elements 2 and 3 is a true value, and is tau₂₃1691.2 us. The simulation result of the static target is shown in fig. 4, wherein fig. 4(a) is the estimation result of the delay difference between the array elements 1 and 2 when the target and the base array are static; fig. 4(b) is the estimation result of delay inequality between array elements 2 and 3 when the target and the base array are static; the open circles in the figure represent the time delay difference without conventional interpolationAnd estimating a simulation result, wherein X represents the delay inequality estimation simulation result obtained by the method. As can be seen from the figure: the conventional non-interpolation method is limited by a sampling rate, the estimation precision of the time delay difference is poor, and the time delay difference is divided into an upper line and a lower line; the method of the invention makes the time delay difference estimation result closer to the true value, and greatly reduces the discrete degree of the time delay difference estimation result. In addition, it can be further appreciated from the figures that: the standard deviation of the delay inequality estimation error of the conventional non-interpolation method is about 49.78us, while the standard deviation of the delay inequality estimation error of the method of the present invention is about 3.728 us. The static target simulation result shows that: when the target and the ternary array are static, the method has better stability.

In order to verify how the application effect of the method is in practical conditions, MATLAB is firstly adopted for numerical simulation analysis, broadband white Gaussian noise is selected as a target radiation signal in the following simulation under the simulation conditions, and the upper and lower limits of the noise frequency are f_l500Hz and f_h1kHz, the array elements are arranged as shown in figure 1, the distance between the array elements is d equal to 50m, the target orientation is theta equal to 60 degrees and 0.1:69.9 degrees, the target distance R is 2000m, the effective sound velocity c is 1500m/s, and the early sampling rate of the system is f_s10kHz, sample length T0.25 s, signal-to-noise ratio snr 20dB, and f is obtained by late 4 times down sampling_s2.5 kHz. Fig. 5 is a simulation result of relative motion between the target and the matrix, where fig. 5(a) is a time delay difference estimation result of array elements 1 and 2 when the target and the matrix move relative to each other; fig. 5(b) is the estimation result of the delay inequality between the array elements 2 and 3 when the target and the matrix move relatively; the hollow circle in the figure represents the time delay difference estimation simulation result under the condition of no conventional interpolation, and the X represents the time delay difference estimation simulation result obtained by the method. As can be seen from fig. 5: the conventional non-interpolation time delay difference estimation result can reflect the movement of a target to a certain extent, but the repeated jump of the estimation value can cause the jump of the ternary array positioning result, and the method can well solve the problem. The moving target simulation results shown in fig. 5 show that: when the target moves at a low speed, the time delay difference estimation method has good continuity.

In conclusion, the high-precision time delay difference estimation method carries out theoretical analysis and experimental verification, combines the advantages of a cross-spectrum method and related peak refinement, and improves the combination process of the cross-spectrum method and the related peak refinement to a certain extent, and adopts related data to directly calculate a time delay difference period value and a time delay difference correction value in a frequency domain without converting the related peak refinement method to a time domain to obtain the time delay period difference value and the time delay difference correction value, so that the calculation amount is reduced, the precision is ensured, and meanwhile, the useful information is ensured not to be lost; the process of the invention is described at 10⁴The sampling rate of Hz level can make the estimation precision of the time delay difference reach 10^-7And s. MATLAB numerical simulation herein has verified the effectiveness of this method.

Compared with a quadratic spline interpolation method, the method has the advantages that the limitation of the quadratic spline interpolation method on local data refinement near a related function peak point is overcome, useful information of sampled data is fully utilized, and the time delay difference estimation precision is improved; compared with a correlation peak refining method, the method avoids the defect that the time delay difference estimation precision is sharply reduced due to cross-correlation operation when the signal-to-noise ratio is too low, can directly calculate the time delay difference period value and the time delay difference correction value in a frequency domain, greatly reduces the operation amount, is convenient for engineering realization, and has strong practicability. However, since the method of the present invention requires phase information, which requires that each array element and the filter amplifier must have better phase consistency, in practical application, the method needs to calibrate the phase of each channel in the ternary array positioning system, and compensate the calibration result to the time delay difference solution of the method of the present invention

Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present invention and are not limited. Although the present invention has been described in detail with reference to the embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (4)

1. A delay inequality estimation method based on phase compensation is used for estimating the delay inequality of different receiving array elements receiving target radiation signals in an underwater acoustic ternary array passive positioning system, and comprises the following steps:

performing phase compensation on the cross-power spectrum signals among the receiving array elements on a frequency domain;

then, energy calculation and peak value search are carried out on the cross-power spectrum signal after phase compensation, and a time delay difference period value is obtained;

then, a time delay difference correction value is obtained for the cross-spectrum signal after the phase compensation according to a least square method at the peak point of the phase compensation;

and finally, adding the delay difference period value and the delay difference correction value to obtain a final delay difference estimation value.

2. The method of phase compensation based delay inequality estimation according to claim 1, further comprising:

step 1), receiving a target radiation signal by each array element in a receiving array, and then preprocessing the target radiation signal to obtain a frequency domain signal of the target radiation signal received by each array element;

step 2), selecting any two receiving array elements, and multiplying frequency domain signals of target radiation signals received by the two receiving array elements after Fourier transformation to obtain frequency domain cross spectrum data of the two receiving array elements;

step 3), performing phase compensation on the frequency domain cross-spectrum data obtained in the step 2);

step 4), energy is obtained for the result obtained in the step 3) after phase compensation, then energy peak value search is carried out according to the result obtained by energy peak value search, and the time value corresponding to the obtained energy peak value point is a delay difference period value;

step 5), solving a phase difference of the frequency domain cross spectrum data after the phase compensation at the energy peak point obtained after the phase compensation, and then solving a time delay difference correction value of the phase difference according to a least square method;

and 6) adding the period value of the time delay difference obtained in the step 4) and the time delay difference correction value obtained in the step 5) to obtain a final time delay difference estimation value of the two receiving array elements selected in the step 2).

3. The method for estimating delay spread based on phase compensation according to claim 2, wherein in the step 1), the preprocessing comprises:

step 1-1), filtering and amplifying target radiation signals received by each array element;

step 1-2), carrying out discrete sampling on the filtered signal at a certain sampling rate;

and 1-3) carrying out Fourier transform on the discrete sampling result to obtain a corresponding frequency domain signal.

4. The method for estimating delay inequality based on phase compensation according to claim 2, wherein in the step 3), the calculation formula for performing phase compensation on the frequency domain cross spectrum data obtained in the step 2) is as follows:

${G^{\&prime;}}_{x_1{x}_2}(\mathrm{\&Delta;}t,f)=G_{x_1{x}_2}\left(f\right)e^{j2\mathrm{\&pi;}\mathrm{\&Delta;}tf},f_l\&le;f\&le;f_h,-\frac{d}{c}\&le;\mathrm{\&Delta;}t=\frac{1}{{\mathrm{kf}}_s}\&le;\frac{d}{c}$

wherein,shows the result of phase compensation of the frequency domain cross-spectrum data,representing frequency domain cross-spectrum data, d is array element spacing, c is effective sound velocity, k is sampling rate lifting multiple during phase compensation, and f_lLower limit of filter, f_hIs the upper limit of the filter, f_sIs the sampling rate.