一种用于水声被动定位的信号到达时延差跟踪算法

摘要

本发明公开了一种用于水声被动定位的信号到达时延差跟踪算法，能够有效地解决基于到达时延差的水声被动定位算法中的多路径干扰问题，且不需要人工指定跟踪算法的到达时延差初始值。包括如下步骤：（1）变量初始化，对信号进行分帧；（2）采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函数值序列；（3）确定当前信号帧到达两个水听器的时延差；（4）将最新的到达时延差估计值作为跟踪初始值；（5）估计新一帧信号的互相关函数值序列，在跟踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值，直至所有信号帧处理完成。

说明书

技术领域

本发明涉及水声信号处理技术领域中的水声被动定位方法，尤其涉及一种用 于水声被动定位的信号到达时延差跟踪算法。

背景技术

水声被动定位是水声信号处理领域的一项重要课题，也是一项难题。目前， 水声被动定位的方法主要有三类：基于到达时延差（TDOA）的方法、目标运动 分析（TMA）方法和匹配场处理方法（MFP）。

基于TDOA的方法是最常用的水声被动定位算法之一，并广泛应用于实时 水声被动定位系统中。该算法根据到达水听器阵列不同基元的信号相对时延来进 行目标方位、距离和深度参数的估计。根据所使用的参数估计方法不同，基于 TDOA的目标参数估计可以分为最大似然估计(MLE)和最小二乘方法(LSE)估计 ^[1]。最小二乘估计不需要任何先验知识，并能给出闭合形式（closed form）的解， 通常用来求解超定方程组，因此在需要实时处理的被动定位系统中得到了广泛的 应用。在无干扰噪声的情况下，根据r维空间中r+1个基元间的r个信号到达时 延差就能够唯一确定目标声源的位置^[2]。有噪声存在时，信号到达时延的测量会 产生误差，为了弥补干扰噪声造成的影响，定位所使用的基元个数要多于r+1个。

在水声被动定位实时处理系统中，TDOA通常首先采用估计两个基元信号的 互相关函数，然后寻找互相关函数的最大值对应时延作为信号的到达时延差^[3]。 Knapp和Cater对互相关法进行了改进，提出了广义互相关法^[4]，以克服互相关 函数估计可能存在的偏差的影响。广义互相关法虽然对TDOA的估计性能有所 改进，但该方法也增加了TDOA估计的计算量。此外，在接收信号中存在经过 多条路径传播到达的信号相叠加时，广义互相关法的性能也受到影响。多路径干 扰导致相关函数中可能出现多个极大值（称为相关峰），由于相关函数的最大值 有可能对应于“伪”相关峰（非同类性质的多径信号相关导致），因此直接通过 寻找互相关函数的最大值计算得到的TDOA通常会发生错误，这种情况在水声 信号传播过程中是普遍存在的。因此，必须采取措施寻找真实的相关峰对应的 TDOA。

参考文献

[1]Huang Y.,Benesty J.,Elko G.W.,Mersereau R.M.Real-Time Passive Source  Localization:A Practical Linear-Correction Least-Squares Approach[J].IEEE  Transaction on Speech and Audio Processing,2001,9(8):943-955.

[2]Smith J.O.,Abel J.S.The Spherical Interpolation Method of Source  Localization[J].IEEE Journal of Oceanic Engineering,1987,12(1):246-252.

[3]张贤达.现代信号处理（第二版）[M].北京：清华大学出版社，2002： 328-330.

[4]Knapp C.H.,Carter G.C.The Generalized Correlation Method for Estimation  of Time Delay[J].IEEE Transaction on Acoustics,Speech and Signal Processing, 1976,24(4):320-327.

发明内容

发明目的：针对现有技术中存在的问题与不足，本发明提供一种用于水声被 动定位的信号到达时延差跟踪算法。在基于TDOA的水声被动定位算法处理过 程中，利用相关函数全局最大值和最小值之间的幅度、时延等相对关系，辨别相 关峰的真伪，从而找到真实的信号到达时延差，并对该时延差进行跟踪，输出稳 定的信号到达时延差估计值。

技术方案：一种用于水声被动定位的信号到达时延差跟踪算法，步骤如下：

（1）变量初始化，对信号进行分帧；

（2）采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函 数值序列；

（3）搜索互相关函数值序列的全局最大值，利用互相关函数值序列的下标 与到达时延差之间的对应关系，确定当前信号帧到达两个水听器的时延差；

（4）重复步骤（2）和（3），直至得到连续5个稳定的到达时延差估计值， 将最新的到达时延差估计值作为跟踪初始值；

（5）使用离散傅里叶变换方法估计新一帧信号的互相关函数值序列，在跟 踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大 值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值；

（6）重复步骤5，直至所有信号帧处理完成。

详细实施步骤为：

1、步骤（1）具体为：

1）定义一个数组TraceArray含有5个元素，初值都置为0；定义计数器变 量Counter，置为0；定义计数器变量FrameCounter，置为1；

2）设两个水听器之间的距离为d（m），水中声速为c（m×s^-1），对两个水 听器接收到的水声信号的采样频率为f_s（Hz），信号长度为每个通道L点，g_L为极大值搜索门限，每帧数据长度为N点，N>0.1f_s且N>(d+g_L)f_s/c+1，帧数M 为不大于L/N的最大正整数；令g_W=d·f_s/c，P为不大于g_W的最大正整数；g_P为 相关函数序列全局最大值与最小值时间差的绝对值门限；

3）数字滤波器H的阶数为k，下限频率为f_L（Hz），上限频率为f_H（Hz）， 且有f_L<f_H<f_s/4，数字滤波器采用FIR滤波器的窗函数法设计。

2、步骤（2）具体为：

1）设各水听器接收到的第FrameCounter帧信号为S_ij，i表示通道序号， i=1,2，j表示数据帧号，j=FrameCounter；

2）对S_ij通过数字滤波器H，得到经过滤波的信号S’_ij；

3）取S’_1j和S’_2j这两个通道的信号，采用互相关函数的傅里叶变换估计法 计算S’_1j和S’_2j两个通道的互相关函数值序列R₂₁(m)，m为整数，且有-(N-1)≤m ≤N-1。

3、步骤（3）具体为：

1）求取R₂₁(m)在下标区间[-(N-1)+P，N-1-P]内的最大值R_MAX，R_MAX对应 的下标为m_MAX，-(N-1)+P≤m_MAX≤N-1-P，对应的信号到达两个水听器的时延 差为m_MAX·c/f_s；

2）求取R₂₁(m)在下标区间[-(N-1)+P，N-1-P]内的最小值R_MIN，R_MIN对应的 下标为m_MIN，-(N-1)+P≤m_MIN≤N-1-P。

4、步骤（4）具体为：

1）如果|m_MAX-m_MIN|≥g_P·f_s，则Counter值加1，并且令TraceArray(Counter)= m_MAX·c/f_s；

2）如果|m_MAX-m_MIN|<g_P·f_s，若|R_MAX|≥|R_MIN|，则Counter值加1，并且令 TraceArray的第Counter个元素TraceArray(Counter)=m_MAXc/f_s；若|R_MAX|<|R_MIN|， 则Counter值置为0；

3）如果Counter值为5，则进行第5）步；如果Counter<5，则FrameCounter 加1，取第FrameCounter帧数据S_ij（i=1,2，j=FrameCounter）；

4）重复步骤（2）和（3），以及步骤（4）的第1）～3）步，直至Counter 值为5；

5）令δ₁=|TraceArray(2)-TraceArray(1)|，δ₂=|TraceArray(3)-TraceArray(2)|， δ₃=|TraceArray(4)-TraceArray(3)|，δ₄=|TraceArray(5)-TraceArray(4)|；

6）若δ₁、δ₂、δ₃、δ₄均小于g_L，Counter置为0，令FrameCounter加1，令 跟踪初始时延差为τ₀=TraceArray(5)/c，进行步骤（5）；

7）若δ₁、δ₂、δ₃、δ₄不均小于g_L，Counter置为0，令FrameCounter加1， 返回至步骤（2）。

5、步骤（5）具体为：

1）设各水听器接收到的第FrameCounter帧信号为S_ij，i表示通道序号， i=1,2，j表示数据帧号，j=FrameCounter；

2）对S_ij通过数字滤波器H，得到经过滤波的信号S’_ij；

3）取S’_1j和S’_2j这两个通道的信号，采用互相关函数的傅里叶变换估计法 计算S’_1j和S’_2j两个通道的互相关函数值序列R₂₁(m)，m为整数，且有-(N-1)≤m ≤N-1；

4）令P₀为不大于τ₀·f_s的最大正整数，Q₀为不大于g_L·f_s/c的最大正整数， 求取R₂₁(m)在下标区间[P₀-Q₀，P₀+Q₀]内的最大值R’_MAX及下标m’_MAX，P₀-Q₀≤m’_MAX≤P₀+Q₀，输出当前帧时延差估计结果τ=m’_MAX/f_s，令τ₀=τ；

6、步骤（6）具体为：

FrameCounter加1，若FrameCounter≤M，则返回至步骤（5）；否则，结束 处理。

为了进一步说明本发明，现将本发明的原理介绍如下：

在水声传播环境中，声源发出的信号到达两个水听器所经过的路径是不同 的，传播所经过的时间延迟也是不同的。当两个水听器距离在一定的范围内，这 两个水听器接收到的信号能够保持较好的相关性，通过寻求两个水听器接收到的 信号的最大相似性发生的时间差，就可以得到信号到达两个水听器的时延差，这 就是基于互相关函数的TDOA估计算法的基本原理。但是，声源发出的信号到 达某个水听器的路径也不是单一的，有可能存在直达路径（D，即声波在传播中 没有和海底或海面发生碰撞）、海面反射路径（S₁，声波与海面碰撞一次后，经 过反射到达水听器）、海底反射路径（B₁，声波与海底碰撞一次后，经过反射到 达水听器）、海面-海底反射路径（S₁B₁，声波先与海面碰撞一次，再与海底碰撞 一次，然后到达水听器）、海底-海面反射路径（B₁S₁，声波先与海底碰撞一次， 再与海面碰撞一次，然后到达水听器）以及多次反射路径（S_nB_n，声波与海底、 海面发生多次碰撞后，到达水听器）等多种传播路径，通常情况下D路径的传 播时延是最小的。经过这些路径传播的信号之间也具有一定的相关性，因此在采 用互相关法估计TDOA时，可能这样的情况：A水听器接收到的信号中的D路 径与B水听器接收信号中的B₁路径相关性较好，导致互相关函数的最大值对应 于B的B₁路径相对于A的D路径的时延差，而真实的传播时延应该是A水听 器接收到的信号中的D路径与B水听器接收信号中的D路径的时延差，这时 TDOA估计就发生了错误。理论分析表明，同种性质的路径相关得到的时延差几 乎相同，其估计误差可以忽略（例如A水听器接收到的信号中的B₁路径与B水 听器接收信号中的B₁路径相关得到的时延差，与两个水听器的D路径相关得到 的时延差可认为是相同的）。

由于海底和海面的界面并不是平整的，经过海底反射路径和海面反射路径传 播的信号幅度并不稳定，因此这种由于不同性质的路径相关得到的到达时延差通 常不是稳定存在的，某些时刻比较明显，某些时刻不明显；而直达路径不和界面 发生作用，经过该路径传播的信号相对比较稳定，直达路径之间相关得到的时延 差是比较稳定的。鉴于水声信号传播的这种特性，可以比较连续几次的TDOA 估计结果，如果这几次的估计值比较接近，则可认为是同种性质的路径相关得到 的，是正确的TDOA估计值。此外，声波经过海面反射后，其相位通常会发生 180度的相移，因此D路径与S₁路径相关得到的值一般小于零，即对应于相关 函数的最小值，根据这一点可以进行D路径与S₁路径相关峰的分辨。

有益效果：与现有技术相比，本发明提供的用于水声被动定位的信号到达时 延差跟踪算法，能够有效地解决基于到达时延差的水声被动定位算法中的多路径 干扰问题，且不需要人工指定跟踪算法的到达时延差初始值，输出的TDOA估 计值稳定。

附图说明

图1为本发明实施例的方法流程图；

图2为水听器A接收的水声信号；

图3为水听器B接收的水声信号；

图4为经过带通滤波的水听器A接收信号；

图5为经过带通滤波的水听器B接收信号；

图6为经过带通滤波的第1帧水听器A接收信号和水听器B接收信号互相 关序列，

图7为经过带通滤波的第2帧水听器A接收信号和水听器B接收信号互相 关序列；

图8为经过跟踪处理的信号到达时延差；

图9为未经过跟踪处理的信号到达时延差。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本 发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发 明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

根据图1，介绍本实施例的用于水声被动定位的信号到达时延差跟踪算法， 包括如下步骤：

（1）变量初始化，对信号进行分帧；

具体为：

1）定义一个数组TraceArray含有5个元素，初值都置为0；定义计数器变 量Counter，置为0；定义计数器变量FrameCounter，置为1；

2）设水听器A和水听器B之间的距离d=22.7（m），水中声速c=1530（m×s^-1）， 对两个水听器接收到的水声信号的采样频率f_s=200（kHz），信号长度为每个通 道L点，L=19064509点，极大值搜索门限g_L=0.3，每帧数据长度为N点，N=40000 点，帧数M=476；g_W=2967.3，P=2967，g_P=0.00025（s）；

3）数字滤波器H的阶数k=512，下限频率f_L=4（kHz），上限频率f_H=12（kHz）； 数字滤波器采用FIR滤波器的窗函数法设计。

（2）采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函 数值序列；水听器A的接收信号用S_1j表示，S₁₁的波形如图2所示，水听器B 的接收信号用S_2j表示，S₂₁的的波形如图3所示；S’₁₁的信号波形如图4所示，S’₂₁的信号波形如图5所示；S’₁₁和S’₂₁的互相关序列R₂₁如图6所示，-39999≤m ≤39999；

（3）搜索互相关函数值序列的全局最大值，利用互相关函数值序列的下标 与到达时延差之间的对应关系，确定当前信号帧到达两个水听器的时延差； R_MAX=475.1，m_MAX=1842，R_MIN=-484.6，m_MIN=275；

（4）重复步骤（2）和（3），直至得到连续5个稳定的到达时延差估计值， 将最新的到达时延差估计值作为跟踪初始值；

|m_MAX-m_MIN|=1570>g_P·f_s=50，则Counter值加1，TraceArray(1)=m_MAX·c/f_s= 14.0913；当Counter值为5时，TraceArray(2)=-14.0925，TraceArray(3)=-14.0886， TraceArray(4)=-14.0896，TraceArray(5)=-14.0949，δ₁=0.0012，δ₂=0.0039， δ₃=0.001，δ₄=0.0053；δ₁、δ₂、δ₃、δ₄均小于g_L，则Counter置为0，FrameCounter 加1，令跟踪初始时延差为τ₀=-0.009212。此处到达时延差为负值，表示信号到 达B水听器的时间早于到达A水听器的时间；

（5）使用离散傅里叶变换方法估计新一帧信号的互相关函数值序列，在跟 踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大 值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值；

S’₁₂和S’₂₂的互相关函数值序列R₂₁如图7所示，P₀=1842，Q₀=39， R’_MAX=524.4，m’_MAX=1842，R’_MIN=-460.0，m’_MIN=1826,τ=m’_MAX/f_s=-0.009212， 令τ₀=τ；

（6）重复步骤5，直至所有信号帧处理完成；

FrameCounter加1，若FrameCounter≤476，则返回至步骤（5）；否则，结 束处理。最后得到的信号到达时延差序列如图8所示，此处为了显示跟踪处理的 效果，将未经跟踪处理的信号到达时延差估计结果也一并画出，如图9所示。