CN104049246A - 一种频率未知的时延差估计方法 - Google Patents

Abstract

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

Description

一种频率未知的时延差估计方法

技术领域

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

背景技术

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

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

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

发明内容

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

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

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

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

每一个频率单元对应的平均时域相关谱乘以相应的加权因子，累加得到最终的 时域相关谱，进而求取最终时延差。 

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

步骤101）对各接收阵元接收的时域目标辐射信号进行傅里叶变换得到相对应的频域信号，其中：m为阵元号，M为阵元数目； 

将任意两个接收阵元接收的目标辐射信号的傅里叶变换后的频域信号X₁(f)和 相乘得到这两个接收阵元的频域互谱数据RX₁₂(f)，采用引导信号与频域互谱数据RX₁₂(f)相乘得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i)，其中：i为频率单元编号且1≤i≤K，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，对接收阵元1和接收阵元2的频域互谱数据RX₁₂(f)中每一个频率单元所对应的时域相关谱R(f_i,τ_j)求取最大值所对应的时间点作为该频率单元的初步时延差估计结果，该初步时延差估计结果记为τ(f_i)，其中，i＝1,...,K，j＝1,...,L，i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数； 

步骤103）更新接收阵元1和接收阵元2接收的时域目标辐射信号，重复进行步骤101)和步骤102)，当重复次数达到预定值N时，得到从起始频率f₁到终止频率f_K，对接收阵元1和接收阵元2的互谱数据RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)，i＝1,...,K，n＝1,...,N，其中：i为频率单元编号，K为终止频率所对应的编号；N取值范围为：5≤N≤15； 

步骤104）计算从起始频率f₁到终止频率f_K中每一个频率单元的时域相关谱R_n(f_i,τ_j)的均值，得到相应的平均时域相关谱；并计算从起始频率f₁到终止频率f_K中 每一个频率单元的N个初步时延差估计结果的方差值，对应结果记为δ_τ(f_i)； 

步骤105）将得到的每一个频率单元的平均相关谱进行时延差方差加权，然后对加权后的各频率单元的相关谱进行累加求和，得到最终的时域相关谱；最后，从最终的时域相关谱中提取出各接收阵元间接收目标辐射信号的时延差，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计。 

上述步骤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}1\≤f\≤f_s.\right.$

其中，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)做逆傅里叶变换，得到从起始频率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)},f_1\≤f_i\≤f_K.\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.$

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

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

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

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

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

步骤104-1）依据步骤103）对重复N次步骤101）得到了从起始频率f₁到终止频率f_K中频率f_i单元的N次时域相关谱R_n(f_i,τ_j)，按下式对频率为f_i的频率单元的N次时域相关谱R_n(f_i,τ_j)取平均，得到频率为f_i的频率单元N次平均时域相关谱，其中，j＝1,...,L，n＝1,...,N； 

$\stackrel{\‾}{R}(f_i,{\mathrm{\τ}}_j)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{R}_n(f_i,{\mathrm{\τ}}_j)}{N}j=1,...,L,n=1,...,N.$

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

$\left\{\begin{array}{c}\mathrm{E\τ}\left(f_i\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\τ}\left(f_i\right)}{N},\\ {\mathrm{\δ}}_{\mathrm{\τ}}\left(f_i\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(\mathrm{\τ}\left(f_i\right)-\mathrm{E\τ}\left(f_i\right))}^2}{N}.\end{array}\right.$

上述步骤105）采用如下公式计算最终的相关谱输出： 

$R_{\mathrm{out}}\left({\mathrm{\τ}}_j\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}\frac{\stackrel{\‾}{R}(f_i,{\mathrm{\τ}}_j)}{{\mathrm{\δ}}_{\mathrm{\τ}}\left(f_i\right)}.j=1,...,L,n=1,...,N.$

其中：R_out(τ_j)为所求的最终的时域相关谱，为所述步骤104-1）所求取的频率为f_i的频率单元的N次平均时域相关谱，δ_τ(f_i)为每个频率单元所对应的时延差方差；i为频率单元编号，K为终止频率所对应的编号，j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数，N取值范围为：5≤N≤15。 

与现有技术相比，本发明的优点在于：提高了时延差估计算法对信噪比的宽容性，即本发明在较低信噪比下可实现水声定位系统中任意两个接收阵元1和2接收目标辐射信号的时延差。 

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

附图说明

图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}1\≤f\≤f_s.\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)},f_1\≤f_i\≤f_K.\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\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)和步骤2)，当重复次数达到预定值N时，可得到从起始频率f₁到终止频率f_K，对RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)，i＝1,...,K，n＝1,...,N。其中：i为频率单元编号，K为终止频率所对应的编号；N根据实际需要而定，一般范围为：5≤N≤15； 

步骤4）、对步骤3）中N次时域相关谱进行平均计算，结果记为j＝1,...,L；其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数。 

步骤5）、对步骤3）中每一个频率单元的时延差估计结果进行方差计算，对应结果记为δ_τ(f_i)，i＝1,...,K，；其中：i为频率单元编号，K为终止频率所对应的编号； 

步骤6）、对步骤4）对每一个频率单元的平均相关谱进行时延差方差加权，然后对加权后的各频率单元的相关谱进行累加求和，如式(1)所示作为最终的相关谱输出。 

$R_{\mathrm{out}}\left({\mathrm{\τ}}_j\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}\frac{\stackrel{\‾}{R}(f_i,{\mathrm{\τ}}_j)}{{\mathrm{\δ}}_{\mathrm{\τ}}\left(f_i\right)}.j=1,...,L,n=1,...,N.---\left(1\right)$

所述步骤1）中，采用引导信号对每一个频率单元求取R(f_i,τ_j)，i＝1,...,K，j＝1,...,L； 

所述步骤2）中，对每一个频率单元相关谱R(f_i,τ_j)，i＝1,...,K，j＝1,...,L输出求取最大值，即每一个频率单元的时延差估计结果，记为τ(f_i)； 

所述步骤4）中，对N次相关谱进行平均计算，结果记为j＝1,...,L； 

所述的步骤5）中，对每一个频率单元的时延差估计结果进行方差计算，对应结果记为δ_τ(f_i)； 

所述的步骤6）中，对每一个频率单元的平均相关谱进行时延差方差加权，然后对加权后的各频率单元的相关谱进行累加求和，如式(1)所示作为最终的相关谱输出。 

在对本发明的方法做详细说明前，首先对本发明方法所适用的接收阵加以描述。图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}1\≤f\≤f_s.\right.---\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)},f_1\≤f_i\≤f_K.\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.---\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)和步骤2)，当重复次数达到预定值N时，可得到从起始频率f₁到终止频率f_K，对RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)，i＝1,...,K，n＝1,...,N。其中：i为频率单元编号，K为终止频率所对应的编号；N根据实际需要而定，一般范围为：5≤N≤15； 

步骤4）、对步骤3）中N次时域相关谱进行平均计算，结果记为j＝1,...,L；其中：i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数。 

步骤5）、对步骤3）中每一个频率单元的时延差估计结果进行方差计算，对应结 果记为δ_τ(f_i)，i＝1,...,K，；其中：i为频率单元编号，K为终止频率所对应的编号； 

步骤6）、对步骤4）对每一个频率单元的平均相关谱进行时延差方差加权，然后对加权后的各频率单元的相关谱进行累加求和，如式(1)所示。 

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

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

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

首先对噪声频率单元进行统计，即对个K-1频率单元进行统计， 

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

即对于噪声的相关谱输出，每个预成时延差输出值是相等的，然后进一步将线谱估计方位的结果累加到(5)式表示结果中，有 

进一步简化有 

$R_{\mathrm{out}}\left({\mathrm{\τ}}_j\right)=(K-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分析，采用引导信号对每一个频率单元进行频域互相关。利用噪声对应频率单元互相关谱的最大值随机，基于目标对应频率单元互相关谱最大值基本一致的特点，统计各频率单元的时延差(TDD)估计结果，理论分析和仿真结果表明：该方法具有较好的有效性，对信噪比的宽容性远好于频域互相关法。该方法为弱线谱时延差估计提供了一个参考思路。本发明涉及一种未知频率的时延差估计方法。 

综上所述，本发明涉及一种未知频率的时延差估计方法，该方法用于估计水声定位系统中不同接收阵元接收目标辐射信号的时延差，所述方法为：求取每一个频率单元的时域相关谱以及相应的时延差；目标辐射信号所占频率单元的时延差为真值时延差比较稳定，非目标辐射信号频率单元对应时延差为非真值时延差比较随机；对每一个频率单元对应的时延差进行统计计算得到该频率单元对应的时延差方差，在对时延差方差进行求倒数得到加权因子；然后，将每一个频率单元对应的平均时域相关谱乘以相应的加权因子，并累加得到最终的时域相关谱；此时，最终时域相关谱中目标辐射信号所对应的频率单元的时域相关谱所占比重比非目标辐射信号所对应的频率单元的时域相关谱所占比重大，进而突出了目标辐射信号所占频率单元对应的时域相关谱在最终时域相关谱中的比值；再从最终的时域相关谱中求取最终的时延差估计值。本方法具有较好的有效性，对信噪比的宽容性远好于频域互相关法，可满足实际工程应用需要。 

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

Claims (6)

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

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

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

每一个频率单元对应的平均时域相关谱乘以相应的加权因子，累加得到最终的时域相关谱，进而求取最终时延差。

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

步骤101）对各接收阵元接收的时域目标辐射信号进行傅里叶变换得到相对应的频域信号，其中：m为阵元号，M为阵元数目；将任意两个接收阵元接收的目标辐射信号的傅里叶变换后的频域信号X₁(f)和相乘得到这两个接收阵元的频域互谱数据RX₁₂(f)，采用引导信号与频域互谱数据RX₁₂(f)相乘得到从起始频率f₁到终止频率f_K中RX₁₂(f)每一频率单元所对应的频域互谱数据RX₁₂(f_i)，其中：i为频率单元编号且1≤i≤K，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，对接收阵元1和接收阵元2的频域互谱数据RX₁₂(f)中每一个频率单元所对应的时域相关谱R(f_i,τ_j)求取最大值所对应的时间点作为该频率单元的初步时延差估计结果，该初步时延差估计结果记为τ(f_i)，其中，i＝1,...,K，j＝1,...,L，i为频率单元编号，K为终止频率所对应的编号；j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数；

步骤103）更新接收阵元1和接收阵元2接收的时域目标辐射信号，重复进行步骤101)和步骤102)，当重复次数达到预定值N时，得到从起始频率f₁到终止频率f_K，对接收阵元1和接收阵元2的互谱数据RX₁₂(f)中每一个频率单元所对应的N个初步的时延差估计结果，记为τ_n(f_i)，i＝1,...,K，n＝1,...,N，其中：i为频率单元编号，K为终止频率所对应的编号；N取值范围为：5≤N≤15；

步骤104）计算从起始频率f₁到终止频率f_K中每一个频率单元的时域相关谱R_n(f_i,τ_j)的均值，得到相应的平均时域相关谱；并计算从起始频率f₁到终止频率f_K中每一个频率单元的N个初步时延差估计结果的方差值，对应结果记为δ_τ(f_i)；

步骤105）将得到的每一个频率单元的平均相关谱进行时延差方差加权，然后对加权后的各频率单元的相关谱进行累加求和，得到最终的时域相关谱；最后，从最终的时域相关谱中提取出各接收阵元间接收目标辐射信号的时延差，完成水声定位系统中不同接收阵元接收目标辐射信号的时延差估计。

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}1\≤f\≤f_s.\right.$

其中，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)做逆傅里叶变换，得到从起始频率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)},f_1\≤f_i\≤f_K.\\ R(f_i,{\mathrm{\τ}}_j)=\mathrm{ifft}\left({\mathrm{RX}}_{12}\left(f_i\right)\right),\end{array}\right.$

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

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

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

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

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

步骤104-1）依据步骤103）对重复N次步骤101）得到了从起始频率f₁到终止频率f_K中频率f_i单元的N次时域相关谱R_n(f_i,τ_j)，按下式对频率为f_i的频率单元的N次时域相关谱R_n(f_i,τ_j)取平均，得到频率为f_i的频率单元N次平均时域相关谱，其中，j＝1,...,L，n＝1,...,N；

$\stackrel{\‾}{R}(f_i,{\mathrm{\τ}}_j)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{R}_n(f_i,{\mathrm{\τ}}_j)}{N}j=1,...,L,n=1,...,N.$

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

$\left\{\begin{array}{c}\mathrm{E\τ}\left(f_i\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\τ}\left(f_i\right)}{N},\\ {\mathrm{\δ}}_{\mathrm{\τ}}\left(f_i\right)=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(\mathrm{\τ}\left(f_i\right)-\mathrm{E\τ}\left(f_i\right))}^2}{N}\end{array}\right.$

根据上式可知：目标辐射信号所占频率单元的时延差方差较小，非目标辐射信号所占频率单元的时延差方差较大；再对时延差方差进行求倒数得到加权因子，则目标辐射信号所占频率单元对应的加权因子较大，非目标辐射信号所占频率单元对应的加权因子较小。

6.根据权利要求2所述的基于频率方差加权的时延差估计方法，其特征在于，所述步骤105）采用如下公式计算最终的相关谱输出：

$R_{\mathrm{out}}\left({\mathrm{\τ}}_j\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}\frac{\stackrel{\‾}{R}(f_i,{\mathrm{\τ}}_j)}{{\mathrm{\δ}}_{\mathrm{\τ}}\left(f_i\right)}.j=1,...,L,n=1,...,N.$

其中：R_out(τ_j)为所求的最终的时域相关谱，为所述步骤104-1）所求取的频率为f_i的频率单元的N次平均时域相关谱，δ_τ(f_i)为每个频率单元所对应的时延差方差；i为频率单元编号，K为终止频率所对应的编号，j时延差对应的时间点编号，L为R(f_i,τ_j)在时间轴上的点数，N取值范围为：5≤N≤15。