CN104237871A - 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

A kind of time delay estimation method based on phase compensation

Technical field

The present invention relates to signal processing field, particularly a kind of time delay estimation method based on phase compensation.

Background technology

Utilizing submarine target institute radiated noise to carry out detection and positioning target is a kind of important means that submarine concealment finds and attacks.Current Passive Location mainly contains ternary array, target motion analysis and Matched-field processing 3 kinds, wherein ternary passive positioning method is still the Passive Location that submarine sonar equipment generally adopts, its key is two time delay estimations between 3 submatrix Received signal strength, and reason is that time delay estimation precision directly determines the positioning precision to target.

Because submarine target institute radiated noise is mostly broadband signal, the cross correlation function of signal all has sharp-pointed relevant peaks main lobe usually, and receiving end can be asked for target emanation signal and arrives delay inequality between each array element by carrying out cross-correlation to received signal.In order to break through the restriction of sampling rate, researcher proposes multiple relevant peaks thinning method, mainly be divided into two types: a class directly carries out refinement method in related function peak value local, as cubic spline interpolation, quadratic spline interpolation, cosine function method of interpolation etc., this class methods operand is little, but precision is low; Another kind of is that process in early stage makes relevant peaks refinement, and as MCZT method, ZOOM-FFT method, BCZT method etc., these class methods are no longer carry out interpolation to peak value local, but realize refinement by frequency domain computing, and operand is large, but precision comparison is high.

In order to improve time delay estimation precision between two paths of signals, J.S. the cross-power spectrum of Bei Endate and A.G. Pi Aisuoer proposition calculating two signals carrys out the delay inequality between estimated signal, but there will be phase fuzzy problem in practical application, tolerance is poor, time delay estimation precision is not high on the contrary.

Summary of the invention

The object of the invention is to time delay estimation method of the prior art and there is phase fuzzy problem in actual applications, tolerance is poor, the not high defect of time delay estimation precision, thus a kind of time delay estimation method that tolerance is high, time delay estimation precision is high is provided.

To achieve these goals, the invention provides a kind of time delay estimation method based on phase compensation, for estimating the different delay inequality receiving array element receiving target radiation signal in underwater sound ternary array Passive Positioning System, described method comprises:

Frequency domain carries out phase compensation to the cross-power spectrum signal between described reception array element;

Then carry out energy to the cross-power spectrum signal after phase compensation to ask for and peak value searching, obtain delay inequality periodic quantity;

Next at the peak point place of phase compensation, by least square method, delay inequality modified value is asked for the cross-spectrum signal after phase compensation;

Finally delay inequality periodic quantity is added with delay inequality modified value, obtains final time delay estimation value.

In technique scheme, described method comprises further:

Each array element receiving target radiation signal in step 1), receiving array, then carries out pre-service to target emanation signal, obtains the frequency-region signal of each array element institute receiving target radiation signal;

Step 2), choose any two and receive array elements, the frequency-region signal received these two after the Fourier transform of target emanation signal that array elements receive is multiplied, and obtains the frequency domain cross-spectrum data that these two receive array elements;

Step 3), to step 2) the frequency domain cross-spectrum data that obtain carry out phase compensation;

Result after step 4), the phase compensation that obtains step 3) is carried out energy and is asked for, and the result then asked for according to energy does energy peak search, and the moment value corresponding to energy peak point obtained is delay inequality periodic quantity;

Step 5), the energy peak point place obtained after phase compensation ask for phase differential to the frequency domain cross-spectrum data after phase compensation, then ask for delay inequality modified value to phase differential according to least square method;

The delay inequality modified value that step 6), delay inequality periodic quantity and the step 5) of step 4) being tried to achieve are tried to achieve is added, and obtains step 2) in selected two final time delay estimation value differences receiving array elements.

In technique scheme, in described step 1), described pre-service comprises:

Step 1-1), filter and amplification is done to the target emanation signal received by each array element;

Step 1-2), with certain sampling rate, discrete sampling is carried out to filtered signal;

Step 1-3), Fourier transform is carried out to discrete sampling result, obtain corresponding frequency-region signal.

In technique scheme, in described step 3), to step 2) the frequency domain cross-spectrum data that the obtain computing formula of carrying out phase compensation is:

${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}$

during phase compensation, sampling rate promotes multiple, f _lfor wave filter lower limit, f _hfor the wave filter upper limit, f _sfor sampling rate.

Strengthen Mutual spectrum tolerance in actual applications, improve time delay estimation precision.

The invention has the advantages that:

The method of the invention can obtain very high-precision delay inequality, and has very strong practicality, effectively can improve underwater sound ternary array Passive Positioning System positioning precision.

Accompanying drawing explanation

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

Fig. 2 is the process flow diagram of the inventive method;

Fig. 3 is distinct methods time delay estimation error to standard deviation;

Fig. 4 (a) for target and basic matrix static time, array element 1,2 time delay estimation result;

Fig. 4 (b) for target and basic matrix static time, array element 2,3 time delay estimation result;

When Fig. 5 (a) is for target and basic matrix relative motion, array element 1,2 time delay estimation result;

When Fig. 5 (b) is for target and basic matrix relative motion, array element 2,3 time delay estimation result.

Embodiment

Now the invention will be further described by reference to the accompanying drawings.

Before method of the present invention is elaborated, first to the inventive method the receiving array that is suitable for be described.

Fig. 1 is the schematic diagram of a certain ternary array, and this ternary array comprises three array elements, represents respectively by numeral 1,2,3, and target is from θ direction radiation signal, and the signal of institute's radiation arrives each array element after underwater acoustic channel is propagated.In order to the convenience illustrated, be described hereinafter for 1,2 array elements to method of the present invention, for other array element, method of the present invention is applicable equally.

Based on all receiving arraies as shown in Figure 1, first the present invention directly carries out phase compensation to the cross-power spectrum signal between array element on frequency domain; Then the cross-power spectrum signal after compensating is carried out that energy is asked for, peak value searching, obtain delay inequality periodic quantity; Next at the peak point place of phase compensation, by least square method, delay inequality modified value is asked for cross-spectrum signal after compensation; Finally delay inequality cycle codomain delay inequality modified value is added and obtains final time delay estimation value.

Below with reference to Fig. 2, specifically comprising of the inventive method:

Each array element receiving target radiation signal in step 1), receiving array, then carries out pre-service to target emanation signal, obtains the frequency-region signal of each array element institute receiving target radiation signal.

The pre-service that target emanation signal does is comprised:

Step 1-1), filter and amplification is done to the target emanation signal received by each array element, filtered signal post is designated as x _m(t), m wherein represents array element number;

Step 1-2), with sample rate f _sto filtered signal x _mt () carries out discrete sampling, obtain x _m(nT _s); T wherein _s=1f _s, represent sampling time interval;

Step 1-3), to discrete sampling result x _m(nT _s) carry out Fourier transform, obtain corresponding frequency domain signal X _m(f), f _l≤ f≤f _h; Wherein, f _lfor wave filter lower limit, f _hfor the wave filter upper limit;

Step 2), the frequency domain signal X that will receive after array element 1 and the Fourier transform receiving the target emanation signal that array element 2 receives ₁(f), f _l≤ f≤f _hwith f _l≤ f≤f _hbe multiplied, obtain the frequency domain cross-spectrum data that these two receive array element f _l≤ f≤f _h, wherein ^*() is complex conjugate.

Step 3), to step 2) the frequency domain cross-spectrum data that obtain f _l≤ f≤f _hcarry out phase compensation.

The computing formula of frequency domain cross-spectrum data being carried out to phase compensation 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}$ Expression frequency domain cross-spectrum data do the result after phase compensation, and d is array element distance, and c is effective sound velocity, and when k is phase compensation, sampling rate promotes multiple.

Result after step 4), the phase compensation that obtains step 3) is carried out energy and is asked for, and the result then asked for according to energy does energy peak search, and the moment value corresponding to energy peak point obtained is delay inequality periodic quantity.

In this step, formula that energy asks for is carried out as shown in Equation (2) to the result after phase compensation:

${{\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)$

Energy peak search computing formula as shown in Equation (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 _12Nrepresent delay inequality periodic quantity.

Step 5), the energy peak point place obtained after phase compensation ask for phase differential to the frequency domain cross-spectrum data after phase compensation, then ask for delay inequality modified value to phase differential according to least square method.

In this step, ask for the formula of phase differential as shown in formula (4) below:

Wherein f _l≤ f≤f _hrepresent phase differential.

τ wherein _12Mrepresent delay inequality modified value.

Step 6), the delay inequality periodic quantity τ that step 4) is tried to achieve _12Nthe delay inequality modified value τ tried to achieve with step 5) _12Mbe added, obtain final time delay estimation value difference

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

Be more than the description to the inventive method basic step, below the necessity of time delay estimation method of the present invention be described.

Hypothetical target radiation signal is bandlimited signal is s (t), and signal length is T, and the frequency limits of signal is respectively f _land f _h.Signal arrives receiving end:

$\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 ₁(t) and x ₂t () is respectively array element 1,2 Received signal strength, τ ₁and τ ₂be respectively the propagation delay that target emanation signal arrives array element 1,2, k ₁and k ₂the amplitude of two paths of signals respectively, n ₁(t) and n ₂t () is respectively white Gaussian noise.

First, do Fourier transform to two paths of signals can 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)$

And then the cross-power spectrum that can obtain between two paths of signals is:

$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, () ^*for conjugate transpose.The phase place that can obtain cross-power spectrum signal according to formula (9) is:

Can be obtained by formula (10): phase place comprise the delay inequality τ of array element 1,2 Received signal strength ₁₂, right adopt least square method to carry out matching, can τ be obtained ₁₂optimal estimation be:

If directly adopted as final time delay estimation, when the distance between array element 1,2 is greater than half wavelength, can cause there is phase ambiguity situation, now, cannot directly by measure the periodic quantity (main value) of array element 1,2 phase of received signal difference, so cannot directly by measure the true delay inequality of array element 1,2 Received signal strength; For this problem, phase compensation array element 1,2 Received signal strength cross-power spectrum can be adopted to solve, then carry out that energy is asked for, peak value searching can obtain array element 1,2 Received signal strength delay inequality periodic quantity (main value) to the cross-power spectrum after compensating, its method is as follows:

Available formula (12) represents phase compensation formula,

$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/ (kf _s): T'] be the phase compensation moment, T' is for compensating cut-off time, f _sfor sampling rate, 1/ (kf _s) be make-up time interval, in formula, k measures the high-precision time-delay difference cycle, similar Interpolation across frequency method to improve.

Then right on each time point f _l≤ f≤f _hask for energy, one group of energy value can be obtained ' Pow (t) is searched for, array element 1,2 delay inequality periodic quantity τ can be measured _12N.Owing to being subject to sampling rate restriction when compensating, delay inequality periodic quantity τ will be caused _12Nand there is deviation between delay inequality true value, this deviation definition is delay inequality modified value by we,

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

Now, τ _12Nvalue ensured to calculate Δ τ _12Mtime without phase ambiguity, namely

Then, right adopt least square method to carry out matching by formula (11), delay inequality modified value τ can be obtained _12Moptimal estimation.

Finally delay inequality periodic quantity and delay inequality modified value are added and can obtain time delay estimation value

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

From above derivation formula, ensureing that formula (14) is set up is the key realizing this method, if τ ₁₂with τ _12Nthe limits of error are:

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

According to formula (14), the adequate condition without phase ambiguity is:

2πf _hτ _12E＜π, (17)

Further simplification can obtain:

f _hτ _12E＜12, (18)

Sampling rate general satisfaction Nyquest sampling thheorem in actual applications, then sample rate f _s> 2f _htime, formula (18) can be rewritten as:

f _sτ _12E＜1, (19)

The physical significance of formula (15) is, if gained energy peak point metrical error is less than a sampled point after phase compensation, and calculation delay difference modified value τ _12Mused there is not phase fuzzy problem.

Sample rate f in systems in practice _s=Nf _h, and in phase compensation formula (12), compensation of phase sampling rate used general larger (i.e. k>=10), then formula (19) can be rewritten as:

f _sτ _12E＜Nk2, (20)

Formula (20) is as long as show that the error of delay inequality periodic quantity is less than time span corresponding to Nk2 sampled point, and this method just can use, and there is not phase fuzzy problem.

Below in conjunction with an example, the implementation procedure of the inventive method is described.

First to received signal filter and amplification is carried out, and with sample rate f _ssampling is carried out to the signal after filter and amplification and obtains x ₁(nT _s) and x ₂(nT _s), then Fourier transform process is carried out to the signal after sampling, conveniently carry out computing, at x ₁(nT _s) and x ₂(nT _s) afterbody carries out benefit 0, makes the length of two sections of sequences be N, the N smallest positive integral that should meet the following conditions here: 1) N > 2 (dfs/c)+1, wherein, d is battle array spacing, and c is the velocity of sound; 2) log ₂n is integer.

To x ₁(nT _s) and x ₂(nT _s) do Fast Fourier Transform (FFT) and obtain:

$\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 _sit is corresponding frequency sampling interval.

If then by M ₁to M ₂between frequency corresponding to frequency all at f _lto f _hin; with be respectively and round up and down.

Secondly, the cross-power spectrum asking for two paths of signals by formula (9) is m ₁≤ m≤M ₂, right back-pushed-type (22) carries out phase compensation to cross-power spectrum.

$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)$

Right ask for energy at each time point, then by detection peak, obtain array element 1,2 delay inequality periodic quantity τ _12N.The phase signal of right back-pushed-type (10) and compensation place of formula (11) process peak point, and then obtain delay inequality modified value τ _12M.

Finally by delay inequality periodic quantity τ _12Nwith delay inequality modified value τ _12Mbe added, final time delay estimation value τ can be obtained ₁₂.

In order to illustrate that the method for the invention compares the nearlyer Cramer-Rao lower limit of conventional method, time delay estimation precision is high, provides the Cramer-Rao lower limit of the time delay estimation standard deviation that this model calculates below.

But G.H. that was generally once pointed out, the Cramer-Rao lower limit of the time delay estimation standard deviation that this model calculates met:

$\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, in the T sampling time, r (f) is x ₁(t) and x ₂(t) coherence spectrum, it is 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 _nnf () is respectively the power spectrum of signal and noise; In simple cases, if K _ss(f) and K _nnf the shape of () is all flat, namely

K _ss(f)＝S ₀,K _nn(f)＝N ₀,f∈(f _l,f _h) (25)

Now, can 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)$

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

Adopt the time delay estimation precision of MATLAB to distinct methods to emulate below, select broadband Gaussian white noise as target emanation signal in emulation, noise frequency bound is respectively f _l=500Hz and f _h=1kHz, array element distance is d=50m, and target azimuth is θ=60 °, and target range is R=2000m, effective sound velocity c=1500m/s, and system sampling rate in early stage is f _s=10kHz, sample length is T=0.25s, and after down-sampled by the later stage 2 times, sampling rate becomes f _s=5kHz.

Fig. 3 is distinct methods time delay estimation error to standard deviation, as can be seen from Figure 3:

1), minimum without interpolation correlation method time delay estimation precision, reason is to limit by sampling rate;

2), with compared with interpolation correlation method, quadratic spline interpolation is carried out to relevant peaks envelope local and breaches sampling rate restriction, time delay estimation precision is greatly improved, but because the process asking for peak envelope in early stage lost a lot of information, so it can not reach very high time delay estimation precision;

3), compared with quadratic spline interpolation, ZOOM-FFT method relevant peaks refinement method uses related data to carry out refinement at frequency domain, and then avoid the limitation of only data near peak envelope being carried out to refinement;

4), compared with additive method, the standard deviation of the time delay estimation error of the method for the invention is minimum, and closer to Cramer-Rao lower limit, especially when low signal-to-noise ratio, the time delay estimation precision of the method for the invention does not significantly reduce; Main cause is that this method combines the advantage of Mutual spectrum and relevant peaks refinement, certain improvement has been done to both cohesive process, this method all adopts related data directly to calculate delay inequality periodic quantity and delay inequality modified value at frequency domain, and do not need to forward time domain to ask for delay inequality periodic quantity and delay inequality modified value as relevant peaks refinement method, thus reduce operand, while guarantee precision, ensure that useful information is without loss.Simulation result shows, under the condition selecting above each parameter, (sampling rate is 10 ³hz), the time delay estimation precision of this method can reach 10 ^-6s.

In order to how verify this method effect in a practical situation, connecing down and first adopting MATLAB to carry out Numerical Simulation Analysis, select broadband Gaussian white noise as target emanation signal during simulated conditions emulates as follows, noise frequency bound is respectively f _l=500Hz and f _h=1kHz, array element distance is d=50m, and target azimuth is θ=60 °, target range R=2000m effective sound velocity c=1500m/s, and system sampling rate in early stage is f _s=10kHz, sample length is T=0.25s, and signal to noise ratio (S/N ratio) is that snr=20dB becomes f by later stage 2 times down-sampled post-sampling rate _s=5kHz.According to above-mentioned simulation analysis condition: the delay inequality between array element 1,2 is true value is τ ₁₂=1631.0us, the delay inequality between array element 2,3 is true value is τ ₂₃=1691.2us.Quiet target simulator result as shown in Figure 4, Fig. 4 (a) wherein for target and basic matrix static time, array element 1,2 time delay estimation result; Fig. 4 (b) for target and basic matrix static time, array element 2,3 time delay estimation result; Open circles in figure indicate without conventional without in interpolation situation time time delay estimation simulation result, " X " represents the time delay estimation simulation result that the method for the invention obtains.As we know from the figure: routine limits by sampling rate without method of interpolation, and time delay estimation ratio of precision is poor, point upper and lower two lines; And the inventive method makes time delay estimation result closer to actual value, greatly reduce the dispersion degree of time delay estimation result.In addition, can learn further from figure: the conventional standard deviation without method of interpolation time delay estimation error is about 49.78us, and the standard deviation of the method for the invention time delay estimation error is about 3.728us.Quiet target simulator result shows: when target and ternary basic matrix are static state, the method for the invention has good stability.

In order to how verify this method effect in a practical situation, connecing down and first adopting MATLAB to carry out Numerical Simulation Analysis, select broadband Gaussian white noise as target emanation signal during simulated conditions emulates as follows, noise frequency bound is respectively f _l=500Hz and f _h=1kHz, array element is installed as shown in Figure 1, and array element distance is d=50m, and target azimuth is θ=60 °: 0.1:69.9 °, target range R=2000m effective sound velocity c=1500m/s, and system sampling rate in early stage is f _s=10kHz, sample length is T=0.25s, and signal to noise ratio (S/N ratio) is that snr=20dB becomes f by later stage 4 times down-sampled post-sampling rate _s=2.5kHz.Fig. 5 be target and basic matrix relative motion time simulation result, when Fig. 5 (a) is wherein for target and basic matrix relative motion, array element 1,2 time delay estimation result; When Fig. 5 (b) is for target and basic matrix relative motion, array element 2,3 time delay estimation result; Open circles in figure indicate without conventional without in interpolation situation time time delay estimation simulation result, " X " represents the time delay estimation simulation result that the method for the invention obtains.As can be known from Fig. 5: routine can react target in motion to a certain extent without method of interpolation time delay estimation result, but its estimated value occurs that saltus step repeatedly can cause ternary array positioning result to occur saltus step, and the method for the invention can well address this problem.Moving-target simulation result shown in Fig. 5 shows: in target with when comparatively low velocity is moved, time delay estimation method of the present invention has good continuity.

In a word, high-precision time-delay difference method of estimation of the present invention has carried out theoretical analysis and experimental verification, high-precision time-delay difference method of estimation of the present invention combines the advantage of Mutual spectrum and relevant peaks refinement, certain improvement has been done to both cohesive process, the present invention all adopts related data directly to calculate delay inequality periodic quantity and delay inequality modified value at frequency domain, and do not need to forward time domain to ask for delay inequality periodic quantity and delay inequality modified value as relevant peaks refinement method, thus reduce operand, while guarantee precision, ensure that useful information is without loss; The method of the invention is 10 ⁴the sampling rate of Hz rank can make time delay estimation precision reach 10 ^-7s.In literary composition, MATLAB numerical simulation has demonstrated the validity of the method all.

The present invention, compared with quadratic spline interpolation, has jumped out the limitation of local data's refinement near quadratic spline interpolation pair correlation function peak point, has taken full advantage of the useful information of sampled data, improve time delay estimation precision; Compared with relevant peaks refinement method, this method avoids signal to noise ratio (S/N ratio) too low time, the shortcoming that the time delay estimation precision that array element Received signal strength causes due to computing cross-correlation sharply declines, and can directly obtain delay inequality periodic quantity and delay inequality modified value at frequency domain, greatly reduce operand, be convenient to Project Realization, there is stronger practicality.But because the method for the invention needs to use phase information, this just requires that each array element and filter amplifier bosom friend must have good phase equalization, so in the application of reality, need to calibrate the phase place of each passage in ternary array positioning system, and during delay inequality calibration result being compensated to the method for the invention resolves

It should be noted last that, above embodiment is only in order to illustrate technical scheme of the present invention and unrestricted.Although with reference to embodiment to invention has been detailed description, those of ordinary skill in the art is to be understood that, modify to technical scheme of the present invention or equivalent replacement, do not depart from the spirit and scope of technical solution of the present invention, it all should be encompassed in the middle of right of the present invention.

Claims (4)

1., based on a time delay estimation method for phase compensation, for estimating the different delay inequality receiving array element receiving target radiation signal in underwater sound ternary array Passive Positioning System, described method comprises:

Frequency domain carries out phase compensation to the cross-power spectrum signal between described reception array element;

Then carry out energy to the cross-power spectrum signal after phase compensation to ask for and peak value searching, obtain delay inequality periodic quantity;

Next at the peak point place of phase compensation, by least square method, delay inequality modified value is asked for the cross-spectrum signal after phase compensation;

Finally delay inequality periodic quantity is added with delay inequality modified value, obtains final time delay estimation value.

2. the time delay estimation method based on phase compensation according to claim 1, it is characterized in that, described method comprises further:

Each array element receiving target radiation signal in step 1), receiving array, then carries out pre-service to target emanation signal, obtains the frequency-region signal of each array element institute receiving target radiation signal;

Step 2), choose any two and receive array elements, the frequency-region signal received these two after the Fourier transform of target emanation signal that array elements receive is multiplied, and obtains the frequency domain cross-spectrum data that these two receive array elements;

Step 3), to step 2) the frequency domain cross-spectrum data that obtain carry out phase compensation;

Result after step 4), the phase compensation that obtains step 3) is carried out energy and is asked for, and the result then asked for according to energy does energy peak search, and the moment value corresponding to energy peak point obtained is delay inequality periodic quantity;

Step 5), the energy peak point place obtained after phase compensation ask for phase differential to the frequency domain cross-spectrum data after phase compensation, then ask for delay inequality modified value to phase differential according to least square method;

The delay inequality modified value that step 6), delay inequality periodic quantity and the step 5) of step 4) being tried to achieve are tried to achieve is added, and obtains step 2) in selected two final time delay estimation value differences receiving array elements.

3. the time delay estimation method based on phase compensation according to claim 2, is characterized in that, in described step 1), described pre-service comprises:

Step 1-1), filter and amplification is done to the target emanation signal received by each array element;

Step 1-2), with certain sampling rate, discrete sampling is carried out to filtered signal;

Step 1-3), Fourier transform is carried out to discrete sampling result, obtain corresponding frequency-region signal.

4. the time delay estimation method based on phase compensation according to claim 2, is characterized in that, in described step 3), to step 2) the frequency domain cross-spectrum data that the obtain computing formula of carrying out phase compensation is:

${G^{\&prime;}}_{x_1{x}_2}(\mathrm{\&Delta;t},f)=G_{x_1{x}_2}\left(f\right)e^{j2\mathrm{\&pi;\&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, ${G^{\&prime;}}_{x_1{x}_2}(\mathrm{\&Delta;t},f),f_l\&le;f\&le;f_h,-\frac{d}{c}\&le;\mathrm{\&Delta;t}=\frac{1}{{\mathrm{kf}}_s}\&le;\frac{d}{c}$ Expression frequency domain cross-spectrum data do the result after phase compensation, represent frequency domain cross-spectrum data, d is array element distance, and c is effective sound velocity, and when k is phase compensation, sampling rate promotes multiple, f _lfor wave filter lower limit, f _hfor the wave filter upper limit, f _sfor sampling rate.