CN104049247A - 一种基于频率方差加权的时延差估计方法 - Google Patents

Abstract

本发明提供了一种基于频率方差加权的时延差估计方法，用于估计水声定位系统中不同接收阵元接收的目标辐射信号的时延差，包含：求取若干次每一个频率单元的时域相关谱以及相应的初始时延差；计算得到的每一个初始时延差对应的频率方差，并对频率方差求倒数形成加权因子；累加每一初值时延差对应的加权因子得到最终时域相关谱，从最终时域相关谱中求取最终时延差。本发明采用引导信号对每一个频率单元进行频域互相关，利用非真值时延差对应频率单元随机，真值时延差对应频率单元基本一致的特点，对每一个时延差估计结果进行频率方差统计，突出真值时延差对应线谱信号在时延差估计中所占比重，弱化非真值时延差对应噪声单元在时延差估计中所占比重。

Description

一种基于频率方差加权的时延差估计方法

技术领域

本发明涉及声纳信号处理领域，特别涉及一种基于频率方差加权的时延差估计方法。 

背景技术

由于主动声呐容易暴露，对水下目标进行探测时常采用被动式声呐。随着减振降噪技术的不断提高，目标辐射噪声相比环境噪声在不断地降低，致使声呐设备对其接收信号所能提供的先验知识也在不断地减少。基于宽带能量积分的信号检测方法在被动声呐检测中已不能满足对水下远程目标的探测需求。 

研究表明：水下目标螺旋桨转动会切割水体产生低频信号，其中部分信号会直接以加性形式出现在目标辐射信号中，部分信号则被船体本身的振动调制到较高的频带。在目标辐射信号中线谱谱级通常连续谱谱级要高出10～25dB，这为实现水下目标远程探测提供了一种可能。 

对此，本发明对目标辐射信号频率未知的情况下，介绍了一种基于频率方差加权的时延差估计方法。对接收信号进行FFT分析，采用引导信号对每一个频率单元进行频域互相关。利用非真值时延差对应频率单元随机，真值时延差对应频率单元基本一致的特点，对每一个时延差估计结果进行频率方差统计计算。理论分析和仿真结果表明：本发明具有较好的有效性和对信噪比的宽容性远好于频域互相关法，本发明为弱线谱时延差估计提供了一个参考思路。 

发明内容

本发明的目的在于，为克服低信噪比下在目标辐射信号频率未知的情况下，现有技术不能很好的实现水声定位系统中阵元间接收目标辐射信号的时延差估计问题，本发明提供了一种基于频率方差加权的时延差估计方法。 

为实现上述目的，本发明提供了一种基于频率方差加权的时延差估计方法，该方法用于估计水声定位系统中不同接收阵元接收的目标辐射信号的时延差，所述方法为： 

求取若干次每一个频率单元的时域相关谱以及相应的初始时延差； 

计算得到的每一个初始时延差对应的频率方差，并对频率方差求倒数形成加权 因子； 

累加每一初值时延差对应的加权因子得到最终时域相关谱，从最终时域相关谱中求取最终时延差。 

上述方法具体包含如下步骤： 

步骤101）对各接收阵元接收的目标辐射信号进行傅里叶变换得到各接收阵元接收的目标辐射信号的频域信号； 

将任意两个接收阵元接收的目标辐射信号的频域信号相乘得到两个接收阵元的频域互谱数据RX₁₂(f)； 

采用引导信号与得到的频域互谱数据RX₁₂(f)相乘得到从起始频率f₁到终止频率f_K中频域互谱数据RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i)，其中：i为频率单元编号，K为终止频率所对应的编号；再对RX₁₂(f_i)做逆傅里叶变换，到从起始频率f₁到终止频率f_K中频域互谱数据RX₁₂(f)每一个频率单元所对应的时域相关谱R(f_i,τ_j)，其中：i＝1,···,K，j＝1,···,L，j为时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数； 

步骤102）从起始频率f₁到终止频率f_K，将频域互谱数据RX₁₂(f)中每一个频率单元所对应的时域相关谱R(f_i,τ_j)的最大值所对应的时间点作为该频率单元的初步时延差估计结果，该初步时延差估计结果记为τ(f_i)，i为频率单元编号且i＝1,···,K，K为终止频率所对应的编号；j为时延差对应的时间点编号，L为时域相关谱R(f_i,τ_j)在时间轴上的点数； 

步骤103）更新各接收阵元接收的时域信号，重复进行步骤101)和步骤102)，当重复次数达到预定值N时，从起始频率f₁到终止频率f_K，对RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)；反向推出：每个初步时延差估计结构对应N个频率单元，记为f_n(τ_j)，其中：i＝1,···,K，n＝1,···,N，j＝1,···,L，i为频率单元编号，K为终止频率所对应的编号；为j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；N取值范围为：5≤N≤15； 

步骤104）将得到的每一个初步时延差估计结果进行频率方差计算，对应结果记为δ_f(τ_j)，其中，j＝1,···,L，j为时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数； 

步骤105）依时间轴为顺序对步骤103）中每个初步时延差估计结果所对应步骤104）所求取的频率方差δ_f(τ_j)进行累计计算，累计结果作为最终的时域相关谱；最后，从最终的时域相关谱中提取出各接收阵元间接收目标辐射信号的时延差估计值，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计。 

上述步骤101）进一步包含： 

步骤101-1）对任意两个接收阵元1和2接收的时域的目标辐射信号x₁(t)和x₂(t)进行傅里叶变换得到x₁(t)和x₂(t)相对应的频域信号X₁(f)和X₂(f)；然后，对X₁(f)和相乘，得到接收阵元1和2的频域互谱数据RX₁₂(f)，具体计算公式为： 

$\left\{\begin{array}{c}{X}_1\left(f\right)=\mathrm{fft}\left(x_1\left(t\right)\right),\\ X_2\left(f\right)=\mathrm{fft}\left(x_2\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f\right)=X_1\left(f\right)X_2^{*}\left(f\right),\end{array}\right.1\≤f\≤f_s.$

其中，fft()为傅里叶变换函数，x₁(t)为接收阵元1接收的目标辐射信号，x₂(t)为接收阵元2的接收的目标信号，f_s为提取接收阵元1和2接收目标辐射信号所用的采样率，RX₁₂(f)表示接收阵元1和2频域信号相乘的结果；()^*表示复共轭函数； 

步骤101-2）采用引导信号与RX₁₂(f)相乘得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i),1≤i≤K，其中：i为频率单元编号，K为终止频率所对应的编号；再对RX₁₂(f_i),1≤i≤K做逆傅里叶变换，得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一个频率单元所对应的时域相关谱R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，其中：j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数，具体计算公式为： 

$\left\{\begin{array}{c}s\left(t\right)=\cos \left(2\mathrm{\π}{f}_{i}t\right),\\ S\left(f\right)=\mathrm{fft}\left(s\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f_i\right)=\frac{\left|S\left(f\right)\right|{\mathrm{RX}}_{12}\left(f\right)}{\max \left(\right|S\left(f\right)\left|\right)},\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.f_1\≤f_i\≤f_K.$

其中，s(t)为引导信号，频率f_i为从起始频率f₁到终止频率f_K中的一个频率单元，i为频率单元编号，K为终止频率所对应的编号。 

上述步骤102）进一步包含： 

步骤102-1）依时间轴上的时间点为顺序对步骤101-2）中的R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，进行逐点比较幅值的大小，并记录最大值所在的位置τ(f_i)，i＝1,···,K，其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数； 

步骤102-2）将步骤102-1）中记录的τ(f_i)作为频率为f_i的频率单元的初步时延差估计结果。 

上述步骤104）进一步包含： 

步骤104-1）依据步骤103）对重复N次步骤101）和步骤102）得到了从起始频率f₁到终止频率f_K中RX₁₂(f)的每一个频率单元所对应的N次初步时延差估计结果τ_n(f_i)，i＝1,···,K，n＝1,···,N；其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；然后反向推出：每个时延差对应N个频率单元，记为f_n(τ_j)，j＝1,···,L，n＝1,···,N，其中：L为时域相关谱R(f_i,τ_j)在时间轴的数据长度；N取值范围为：5≤N≤15，对f_n(τ_j)按如下公式求取每个初步时延差估计结果所对应的频率方差： 

$\left\{\begin{array}{c}\mathrm{Ef}\left({\mathrm{\τ}}_j\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{f}_n\left({\mathrm{\τ}}_j\right)}{N},\\ {\mathrm{\δ}}_f\left({\mathrm{\τ}}_j\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(f_n\left({\mathrm{\τ}}_j\right)-\mathrm{Ef}\left({\mathrm{\τ}}_j\right))}^2}{N}.\end{array}\right.$

上述步骤105）进一步包含： 

步骤105-1）依时间轴为顺序对步骤104）求取所得每个时延差对应的频率方差δ_f(τ_j)进行累计计算，即当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则最终所求的时域相关谱R_out(f_n(τ_j))需累加该时延差估计结果对应的频率方差的倒数作为最后的输出值，具体计算公式为： 

R_out(f_n(τ_j))＝R_out(f_n(τ_j))+1/δ_f(τ_j),j＝1,···,L；n＝1,···,N. 

上式具体过程如下：首先将最终所求的时域相关谱R_out(f_n(τ_j))=0，当从第1次开始依次到第N次结束，当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则需将R_out(f_n(τ_j))累加该时延差估计结果对应的频率方差的倒数，第N次结束即得到最终所求的时域相关谱R_out(f_n(τ_j))，进而可从R_out(f_n(τ_j))求出各阵元间接收信号的时延差，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计，其中：L为时域相关谱R(f_i,τ_j)在时间轴的数据长度，N取值范围为：5≤N≤15。 

与现有技术相比，本发明的优点在于：提高了时延差估计算法对信噪比的宽容性，即本发明在较低信噪比下可实现阵元1、2接收目标辐射信号的时延差。 

本方法在目标辐射信号具有稳定线谱的情况下，对接收信号进行FFT分析，采用引导信号对每一个频率单元进行频域互相关。利用非真值时延差对应频率单元随机，真值时延差对应频率单元基本一致的特点，对每一个时延差估计结果进行频率方差统计计算，进而突出真值时延差对应的线谱信号在时延差估计中所占比重，弱化非真值时延差对应的噪声单元在时延差估计中所占比重。本方法具有较好的有效性，对信噪比的宽容性远好于频域互相关法，可满足实际工程应用需要。 

附图说明

图1(a)是本发明提供的基于本发明所述算法的定位系统虚拟装置，虚线包含部分为本发明所述部分； 

图1(b)是本发明提供的基于频率方差加权的时延差估计方法流程图； 

图2本发明的接收信号阵元布阵图； 

图3(a)为目标辐射信号对应的频谱分析结果； 

图3(b)为阵元1接收信号的频谱分析结果； 

图4(a)为频域互相关时延差估计结果； 

图4(b)为频域互相关时延差估计结果局部放大； 

图5(a)为本发明方法的时延差估计结果； 

图5(b)为本发明方法的时延差估计结果局部放大； 

图6不同信噪比下，采用现有技术和采用本发明的方案两种方法的时延差估计概率对比图。 

具体实施方式

下面结合附图对本发明的技术方案进行进一步的描述。 

图1(a)是本发明提供的基于本发明所述算法的定位系统虚拟装置，虚线包含部分为本发明在现有的水下目标定位系统中的所处的位置。 

本发明提供了一种基于频率方差加权的时延差估计方法，包括： 

步骤1）、以阵元1、2为例，对阵元1、2接收的时域信号x₁(t)和x₂(t)进行傅里叶变换得到x₁(t)和x₂(t)相对应的频域信号X₁(f)和X₂(f)；然后，对X₁(f)和相乘，得到阵元1、2的频域互谱数据RX₁₂(f)，具体计算公式为： 

$\left\{\begin{array}{c}{X}_1\left(f\right)=\mathrm{fft}\left(x_1\left(t\right)\right),\\ X_2\left(f\right)=\mathrm{fft}\left(x_2\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f\right)=X_1\left(f\right)X_2^{*}\left(f\right),\end{array}\right.1\≤f\≤f_s.$

其中，fft()为傅里叶变换函数，x₁(t)为阵元1的接收信号，x₂(t)为阵元2的接收信号，f_s为提取阵元1、2接收时域信号的采样率，RX₁₂(f)表示阵元1、2频域信号相乘的结果；()^*表示复共轭函数； 

依据表示阵元1、2频域信号相乘结果，采用引导信号与RX₁₂(f)相乘，可得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i),1≤i≤K，其中：i为频率单元编号，K为终止频率所对应的编号。再对RX₁₂(f_i),1≤i≤K做逆傅里叶变换，可得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一个频率单元所对应的时域相关谱R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，其中：j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数。具体计算公式为： 

$\left\{\begin{array}{c}s\left(t\right)=\cos \left(2\mathrm{\π}{f}_{i}t\right),\\ S\left(f\right)=\mathrm{fft}\left(s\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f_i\right)=\frac{\left|S\left(f\right)\right|{\mathrm{RX}}_{12}\left(f\right)}{\max \left(\right|S\left(f\right)\left|\right)},\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.f_1\≤f_i\≤f_K.$

其中，s(t)为引导信号（需为单频信号，此处用余弦函数表示，时间长度和阵元1、2接收信号一样，初始相位为0，幅度无要求），频率f_i为从起始频率f₁到终止频率f_K中的一个频率单元，i为频率单元编号，K为终止频率所对应的编号 

步骤2）、依时间轴上的时间点为顺序对步骤1）中的R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，进行逐点比较幅值的大小，并记录最大值所在的位置τ(f_i)，i＝1,···,K即可作为本次频率f_i单元的时延差估计结果。其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数。 

步骤3）、更新各接收阵元接收的时域信号x_m(t),m＝1,2,M，重复进行步骤1)和步骤2)，当重复次数达到预定值N时，可得到从起始频率f₁到终止频率f_K，对RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)，i＝1,···,K，n＝1,···,N；反向推出：每个时延差对应N个频率单元，记为f_n(τ_j)，j＝1,···,L，n＝1,···,N。其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；N根据实际需要而定，一般范围为：5≤N≤15； 

步骤4）将得到的每一个时延差估计结果进行频率方差计算，对应结果记为δ_f(τ_j)，j＝1,···,L，j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数； 

步骤5）依时间轴为顺序对步骤4）求取所得每个时延差对应的频率方差δ_f(τ_j)进行累计计算，即当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则本发明最终所求的时域相关谱R_out(f_n(τ_j))需累加该时延差估计结果对应的频率方差的倒数 作为最后的输出值，具体计算公式为： 

R_out(f_n(τ_j))＝R_out(f_n(τ_j))+1/δ_f(τ_j),j＝1,···,L；n＝1,···,N.      (1) 

上式具体过程如下：首先将本发明最终所求的时域相关谱R_out(f_n(τ_j))=0，当从第1次开始依次到第N次结束，当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则需将R_out(f_n(τ_j))累加该时延差估计结果对应的频率方差的倒数，第N次结束可得到本发明最终所求的时域相关谱R_out(f_n(τ_j))，进而可从R_out(f_n(τ_j))求出各阵元间接收信号的时延差，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计。其中：j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；N根据实际需要而定，一般范围为：5≤N≤15。 

实施例： 

首先对本发明方法所适用的接收阵加以描述，图2为一两接收阵元的示意图，目标从θ方向辐射信号，经水声信道传播后到达阵元。以该两阵元为例，下面对本发明的方法做详细说明。 

参考图2本发明的方法包括以下步骤： 

步骤1）、对阵元1、2接收的时域信号x₁(t)和x₂(t)进行傅里叶变换得到x₁(t)和x₂(t)相对应的频域信号X₁(f)和X₂(f)；然后，对X₁(f)和相乘，得到阵元1、2的频域互谱数据RX₁₂(f)，具体计算公式为： 

$\left\{\begin{array}{c}{X}_1\left(f\right)=\mathrm{fft}\left(x_1\left(t\right)\right),\\ X_2\left(f\right)=\mathrm{fft}\left(x_2\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f\right)=X_1\left(f\right)X_2^{*}\left(f\right),\end{array}\right.1\≤f\≤f_s.---\left(2\right)$

其中，fft()为傅里叶变换函数，x₁(t)为阵元1的接收信号，x₂(t)为阵元2的接收信号，f_s为提取阵元1、2接收时域信号的采样率，RX₁₂(f)表示阵元1、2频域信号相乘的结果；()^*表示复共轭函数； 

依据表示阵元1、2频域信号相乘结果，采用引导信号与RX₁₂(f)相乘，可得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i),1≤i≤K，其中：i为频率单元编号，K为终止频率所对应的编号。再对 RX₁₂(f_i),1≤i≤K做逆傅里叶变换，可得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一个频率单元所对应的时域相关谱R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，其中：j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数。具体计算公式为： 

$\left\{\begin{array}{c}s\left(t\right)=\cos \left(2\mathrm{\π}{f}_{i}t\right),\\ S\left(f\right)=\mathrm{fft}\left(s\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f_i\right)=\frac{\left|S\left(f\right)\right|{\mathrm{RX}}_{12}\left(f\right)}{\max \left(\right|S\left(f\right)\left|\right)},\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.f_1\≤f_i\≤f_K.---\left(3\right)$

其中，s(t)为引导信号（需为单频信号，此处用余弦函数表示，时间长度和阵元1、2接收信号一样，初始相位为0，幅度无要求），频率f_i为从起始频率f₁到终止频率f_K中的一个频率单元，i为频率单元编号，K为终止频率所对应的编号 

步骤2）、依时间轴上的时间点为顺序对步骤1）中的R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，进行逐点比较幅值的大小，并记录最大值所在的位置τ(f_i)，i＝1,···,K即可作为本次频率f_i单元的时延差估计结果。其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数。 

步骤3）、更新各接收阵元接收的时域信号x_m(t),m＝1,2,M，重复进行步骤1)和步骤1)，当重复次数达到预定值N时，可得到从起始频率f₁到终止频率f_K，对RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)，i＝1,···,K，n＝1,···,N；反向推出：每个时延差对应N个频率单元，记为f_n(τ_j)，j＝1,···,L，n＝1,···,N。其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；N根据实际需要而定，一般范围为：5≤N≤15； 

步骤4）将得到的每一个时延差估计结果进行频率方差计算，对应结果记为δ_f(τ_j)，j＝1,···,L，j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数； 

步骤5）依时间轴为顺序对步骤4）求取所得每个时延差对应的频率方差δ_f(τ_j)进行累计计算，即当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则本发明最终所求的时域相关谱R_out(f_n(τ_j))需累加该时延差估计结果对应的频率方差的倒数作为最后的输出值，具体计算公式为如式(1)所示。 

以上是对本发明方法基本步骤的描述，下面对这些步骤做进一步的说明。 

在步骤5）中，对基于频率方差加权实现未知频率时延差估计法。下面对这一基于频率方差加权的必要性进行陈述。 

下面给出相关谱输出中，阵元1、2接收目标信号时延差的输出值与其它时延差输出值的大小关系。设最小和最大预成时延差为τ₁、τ_L，最小和最大频率单元为f₁、f_K，对阵元1、2接收目标信号真值时延差对应的频率单元为为，共进行N帧信号统计。假设每一个时延差估计结果对应的频率单元均服从均匀分布，非真值时延差和真值时延差对应的频率方差分别为δ_n、δ_s，有 

首先对非真值时延差进行统计，即对个L-1时延差估计结果进行统计 

$R_{\mathrm{out}}\left({\mathrm{\τ}}_j\right)=(L-1)\·\frac{N}{f_K-f_1}\·\frac{1}{{\mathrm{\δ}}_n},j=\mathrm{1,2},\·\·\·,L.---\left(5\right)$

即对于非真值时延差的相关谱输出，每个时延差幅值是相等的，然后进一步将真值时延差对应频率单元的方差结果累加到(5)式表示结果中，有 

进一步简化有 

$R_{\mathrm{out}}\left({\mathrm{\τ}}_j\right)=(L-1)\·\frac{N}{\sqrt{12}{\mathrm{\δ}}_n^2}+\frac{N}{{\mathrm{\δ}}_s^2}.---\left(7\right)$

因此当真值时延差对应频率单元方差很小时，即每帧时延差估计结果对应频率 均接近于阵元1、2接收信号的线谱信号频率，即δ_s值较小，比较(5)式和(7)式可以看出阵元1、2接收信号的真值时延差附近的相关谱出值将远远大于其它时延差相关谱输出值。 

仿真分析基于远场平面波进行，假设目标相对于阵元1、2的方位为60°，辐射信号包括高斯带限白噪声和线谱成份，白噪声带宽为60-300Hz，线谱频率为100Hz，线谱谱级与白噪声平均谱级比为18dB。目标辐射信号对应的频谱如图3(a)所示。干扰为带限白噪声，目标辐射噪声谱级与干扰噪声谱级比为-16dB，则阵元1、2接收信号的线谱谱级与干扰噪声平均谱级比为2dB。 

图3(a)为目标辐射信号对应的频谱分析结果； 

图3(b)为阵元1接收信号的频谱分析结果。 

从图3(a)可以明显看出目标在频率100Hz处由强线谱存在，而我们不能从图3(b)中得到目标辐射的线谱信号位置，如果采用常规频域相关法求取时延差，干扰噪声与线谱信号在整个相关谱中的比重一样，即非真值时延差与真值时延差在整个相关谱中的比重一样，低信噪比下不能实现时延差估计。 

现假定阵元1、2间距为15m，有效声速为c＝1500m/s，对先前设定的信号按采样率为f_s＝2500Hz，由此可得阵元1、2接收信号的时延差点数为τ＝d·f_s·cos(π/3)/c≈13；现一次采集长度为T＝10s阵元1、2拾取数据，对采集数据分10段，每段分240个频带进行按2种方法进行了仿真分析：第1种方法是基于频域互相关法，即整个频带采用互相关对阵元间信号时延差求取；第2种方法是基于频率方差加权法，即对每一个时延差估计结果统计对应频率单元方差，再进行基于频率方差加权的时延差统计。仿真结果如图4、图5所示。 

图4(a)为频域互相关时延差估计结果； 

图4(b)为频域互相关时延差估计结果局部放大； 

图5(a)为本发明方法的时延差估计结果； 

图5(b)为本发明方法的时延差估计结果局部放大。 

从图4(b)可以看出，基于频域互相关法不能实现对阵元1、2接收信号时延差估计，而从图5(b)互相谱中可以有效地实现对阵元1、2接收信号时延差估计(由于MATLAB仿真中是从1开始，所以图中14即为13)。此仿真表明：在此条件下，本 发明相比频域互相关法可较好的实现阵元1、2接收信号的时延差估计。 

为了验证不同信噪比下，本发明与频域互相关时延差估计概率如何，采用如下的仿真条件进行MATLAB仿真，目标相对于阵元1、2的方位为60°，辐射信号只有线谱成份，干扰噪声带宽为60-300Hz，线谱频率为100Hz。图6给出了线谱与干扰噪声信噪比为SNR＝-28-10dB时，本发明方法与频域互相关的时延差估计概率。 

图6不同信噪比下，两种方法的时延差估计概率 

从图6中可以得到，本发明方法在SNR＝-16dB时的时延差估计概率大于50%，而频域互相关法在SNR＝-7dB时的时延差估计概率已小于50%；且在SNR＝-16--5dB时，本发明方法相比频域互相关法的时延估计概率高出40%。图6的仿真结果表明本方法时延差估计对信噪比的宽容性远好于频域互相关法。

总之，本发明涉及一种基于频率方差加权的时延差估计方法。研究表明：水下目标螺旋桨转动会切割水体产生低频信号，其中部分信号会直接以加性形式出现在目标辐射信号中，部分信号则被船体本身的振动调制到较高的频带。在目标辐射信号中线谱谱级通常连续谱谱级要高出10～25dB，这为实现水下目标远程定位提供了一种可能。对此，本发明对目标辐射信号频率未知的情况下，介绍了一种基于频率方差加权的时延差估计方法。对接收信号进行FFT分析，采用引导信号对每一个频率单元进行频域互相关。利用非真值时延差对应频率单元随机，真值时延差对应频率单元基本一致的特点，对每一个时延差估计结果进行频率方差统计计算。理论分析和仿真结果表明：本发明具有较好的有效性和对信噪比的宽容性远好于频域互相关法，本发明为弱线谱时延差估计提供了一个有效的技术方案。 

综上所述，本发明提供了一种基于频率方差加权的时延差估计方法，该方法用于估计水声定位系统中不同接收阵元接收的目标辐射信号的时延差，所述方法为：求取每一个频率单元的时域相关谱以及相应的时延差；将目标辐射信号所占频率单元的时延差作为真值时延差，将非目标辐射信号频率单元对应时延差作为非真值时延差，进行如下反向推导：真值时延差对应的频率单元为目标辐射信号频率单元比较稳定，非真值时延差对应频率单元比较随机；即真值时延差对应的频率方差较小，非真值时延差对应的频率方差较大；对每一时延差对应频率单元形成频率方差求倒数形成加权因子，此时真值时延差对应的加权因子较大，非真值时延差对应的加权因子较小；通过将最终所求时域相关谱中的每一初值时延差位置处幅值置为加权因子，无初值时延差位置处幅值置为0；此时最终时域相关谱中真值时延差处的幅值较 大，进而突出了目标辐射信号所占频率单元对应的时延差幅值，以便我们从最终的时域相关谱中求取时延差。本发明在较低信噪比下可实现任意两个接收阵元接收目标辐射信号的时延差。 

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。 

Claims (6)

1.一种基于频率方差加权的时延差估计方法，该方法用于估计水声定位系统中不同接收阵元接收的目标辐射信号的时延差，所述方法为：

求取若干次每一个频率单元的时域相关谱以及相应的初始时延差；

计算得到的每一个初始时延差对应的频率方差，并对频率方差求倒数形成加权因子；

累加每一初值时延差对应的加权因子得到最终时域相关谱，从最终时域相关谱中求取最终时延差。

2.根据权利要求1所述的基于频率方差加权的时延差估计方法，其特征在于，所述方法具体包含如下步骤：

步骤101）对各接收阵元接收的目标辐射信号进行傅里叶变换得到各接收阵元接收的目标辐射信号的频域信号；将任意两个接收阵元接收的目标辐射信号的频域信号相乘得到两个接收阵元的频域互谱数据RX₁₂(f)；

采用引导信号与得到的频域互谱数据RX₁₂(f)相乘得到从起始频率f₁到终止频率f_K中频域互谱数据RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i)，其中：i为频率单元编号，K为终止频率所对应的编号；再对RX₁₂(f_i)做逆傅里叶变换，到从起始频率f₁到终止频率f_K中频域互谱数据RX₁₂(f)每一个频率单元所对应的时域相关谱R(f_i,τ_j)，其中：i＝1,···,K，j＝1,···,L，j为时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；

步骤102）从起始频率f₁到终止频率f_K，将频域互谱数据RX₁₂(f)中每一个频率单元所对应的时域相关谱R(f_i,τ_j)的最大值所对应的时间点作为该频率单元的初步时延差估计结果，该初步时延差估计结果记为τ(f_i)，i为频率单元编号且i＝1,···,K，K为终止频率所对应的编号；j为时延差对应的时间点编号，L为时域相关谱R(f_i,τ_j)在时间轴上的点数；

步骤103）更新各接收阵元接收的时域信号，重复进行步骤101)和步骤102)，当重复次数达到预定值N时，从起始频率f₁到终止频率f_K，对RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)；反向推出：每个初步时延差估计结构对应N个频率单元，记为f_n(τ_j)，其中：i＝1,···,K，n＝1,···,N，j＝1,···,L，i为频率单元编号，K为终止频率所对应的编号；为j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；N取值范围为：5≤N≤15；

步骤104）将得到的每一个初步时延差估计结果进行频率方差计算，对应结果记为δ_f(τ_j)，真值时延差对应的频率方差较小，非真值时延差对应的频率方差较大；其中，j＝1,···,L，j为时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；

步骤105）依时间轴为顺序对步骤103）中每个初步时延差估计结果所对应步骤104）所求取的频率方差δ_f(τ_j)进行求倒数形成加权因子，此时真值时延差对应的加权因子较大，非真值时延差对应的加权因子较小；通过将最终所求时域相关谱中的每一初值时延差位置处幅值置为加权因子，无初值时延差位置处幅值置为0；最后，从最终的时域相关谱中提取出各接收阵元间接收目标辐射信号的时延差估计值，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计。

3.根据权利要求2所述的基于频率方差加权的时延差估计方法，其特征在于，所述步骤101）进一步包含：

步骤101-1）对任意两个接收阵元1和2接收的时域的目标辐射信号x₁(t)和x₂(t)进行傅里叶变换得到x₁(t)和x₂(t)相对应的频域信号X₁(f)和X₂(f)；然后，对X₁(f)和相乘，得到接收阵元1和2的频域互谱数据RX₁₂(f)，具体计算公式为：

$\left\{\begin{array}{c}{X}_1\left(f\right)=\mathrm{fft}\left(x_1\left(t\right)\right),\\ X_2\left(f\right)=\mathrm{fft}\left(x_2\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f\right)=X_1\left(f\right)X_2^{*}\left(f\right),\end{array}\right.1\≤f\≤f_s.$

其中，fft()为傅里叶变换函数，x₁(t)为接收阵元1接收的目标辐射信号，x₂(t)为接收阵元2的接收的目标信号，f_s为提取接收阵元1和2接收目标辐射信号所用的采样率，RX₁₂(f)表示接收阵元1和2频域信号相乘的结果；()^*表示复共轭函数；

步骤101-2）采用引导信号与RX₁₂(f)相乘得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i),1≤i≤K，其中：i为频率单元编号，K为终止频率所对应的编号；再对RX₁₂(f_i),1≤i≤K做逆傅里叶变换，得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一个频率单元所对应的时域相关谱R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，其中：j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数，具体计算公式为：

$\left\{\begin{array}{c}s\left(t\right)=\cos \left(2\mathrm{\π}{f}_{i}t\right),\\ S\left(f\right)=\mathrm{fft}\left(s\left(t\right)\right),\\ {\mathrm{RX}}_{12}\left(f_i\right)=\frac{\left|S\left(f\right)\right|{\mathrm{RX}}_{12}\left(f\right)}{\max \left(\right|S\left(f\right)\left|\right)},\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.f_1\≤f_i\≤f_K.$

其中，s(t)为引导信号，频率f_i为从起始频率f₁到终止频率f_K中的一个频率单元，i为频率单元编号，K为终止频率所对应的编号。

4.根据权利要求2所述的基于频率方差加权的时延差估计方法，其特征在于，所述步骤102）进一步包含：

步骤102-1）依时间轴上的时间点为顺序对步骤101-2）中的R(f_i,τ_j)，i＝1,···,K，j＝1,···,L，进行逐点比较幅值的大小，并记录最大值所在的位置τ(f_i)，i＝1,···,K，其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；

步骤102-2）将步骤102-1）中记录的τ(f_i)作为频率为f_i的频率单元的初步时延差估计结果。

5.根据权利要求2所述的基于频率方差加权的时延差估计方法，其特征在于，所述步骤104）进一步包含：

步骤104-1）依据步骤103）对重复N次步骤101）和步骤102）得到了从起始频率f₁到终止频率f_K中RX₁₂(f)的每一个频率单元所对应的N次初步时延差估计结果τ_n(f_i)，i＝1,···,K，n＝1,···,N；其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；然后反向推出：每个时延差对应N个频率单元，记为f_n(τ_j)，j＝1,···,L，n＝1,···,N，其中：L为时域相关谱R(f_i,τ_j)在时间轴的数据长度；N取值范围为：5≤N≤15，对f_n(τ_j)按如下公式求取每个初步时延差估计结果所对应的频率方差：

$\left\{\begin{array}{c}\mathrm{Ef}\left({\mathrm{\τ}}_j\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{f}_n\left({\mathrm{\τ}}_j\right)}{N},\\ {\mathrm{\δ}}_f\left({\mathrm{\τ}}_j\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(f_n\left({\mathrm{\τ}}_j\right)-\mathrm{Ef}\left({\mathrm{\τ}}_j\right))}^2}{N}.\end{array}\right.$

根据上式可知：在真值时延差对应的频率单元为目标辐射信号频率单元比较稳定，非真值时延差对应频率单元比较随机时，可得真值时延差对应的频率方差较小，非真值时延差对应的频率方差较大。

6.根据权利要求2所述的基于频率方差加权的时延差估计方法，其特征在于，所述步骤105）进一步包含：

步骤105-1）依时间轴为顺序对步骤104）求取所得每个时延差对应的频率方差δ_f(τ_j)进行求倒数形成加权因子，此时真值时延差对应的加权因子较大，非真值时延差对应的加权因子较小；通过将最终所求时域相关谱中的每一初值时延差位置处幅值置为加权因子，无初值时延差位置处幅值置为0，即当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则最终所求的时域相关谱R_out(f_n(τ_j))需累加该时延差估计结果对应的频率方差的倒数作为最后的输出值，具体计算公式为：

R_out(f_n(τ_j))＝R_out(f_n(τ_j))+1/δ_f(τ_j),j＝1,···,L；n＝1,···,N.

上式具体过程如下：首先将最终所求的时域相关谱R_out(f_n(τ_j))=0，当从第1次开始依次到第N次结束，当某一时延差估计结果f_n(τ_j)对应频率方差为δ_f(τ_j)时，则需将R_out(f_n(τ_j))累加该时延差估计结果对应的频率方差的倒数，第N次结束即得到最终所求的时域相关谱R_out(f_n(τ_j))，进而可从R_out(f_n(τ_j))求出各阵元间接收信号的时延差，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计，其中：L为时域相关谱R(f_i,τ_j)在时间轴的数据长度，N取值范围为：5≤N≤15。