专利名称：一种用于水声被动定位的信号到达时延差跟踪算法的制作方法
技术领域：
本发明涉及水声信号处理技术领域中的水声被动定位方法，尤其涉及一种用于水声被动定位的信号到达时延差跟踪算法。
背景技术：
水声被动定位是水声信号处理领域的一项重要课题，也是一项难题。目前，水声被动定位的方法主要有三类:基于到达时延差(TDOA)的方法、目标运动分析(TMA)方法和匹配场处理方法(MFP)。基于TDOA的方法是最常用的水声被动定位算法之一，并广泛应用于实时水声被动定位系统中。该算法根据到达水听器阵列不同基元的信号相对时延来进行目标方位、距离和深度参数的估计。根据 所使用的参数估计方法不同，基于TDOA的目标参数估计可以分为最大似然估计(MLE)和最小二乘方法(LSE)估计[1]。最小二乘估计不需要任何先验知识，并能给出闭合形式(closed form)的解，通常用来求解超定方程组，因此在需要实时处理的被动定位系统中得到了广泛的应用。在无干扰噪声的情况下，根据r维空间中r+Ι个基元间的r个信号到达时延差就能够唯一确定目标声源的位置[2]。有噪声存在时，信号到达时延的测量会产生误差，为了弥补干扰噪声造成的影响，定位所使用的基元个数要多于r+Ι个。在水声被动定位实时处理系统中，TDOA通常首先采用估计两个基元信号的互相关函数，然后寻找互相关函数的最大值对应时延作为信号的到达时延差 ]。Knapp和Cater对互相关法进行了改进，提出了广义互相关法[4]，以克服互相关函数估计可能存在的偏差的影响。广义互相关法虽然对TDOA的估计性能有所改进，但该方法也增加了 TDOA估计的计算量。此外，在接收信号中存在经过多条路径传播到达的信号相叠加时，广义互相关法的性能也受到影响。多路径干扰导致相关函数中可能出现多个极大值(称为相关峰)，由于相关函数的最大值有可能对应于“伪”相关峰(非同类性质的多径信号相关导致)，因此直接通过寻找互相关函数的最大值计算得到的TDOA通常会发生错误，这种情况在水声信号传播过程中是普遍存在的。因此，必须采取措施寻找真实的相关峰对应的TD0A。参考文献[I]Huang Y., Benesty J., Elko G.ff., Mersereau R.M.Real-Time PassiveSource Localization:A Practical Linear-Correction Least-Squares Approach[J].1EEE Transaction on Speech and Audio Processing, 2001, 9(8):943-955.
[2]Smith J.0., Abel J.S.The Spherical Interpolation Method of SourceLocalization[J].1EEE 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 forEstimation of Time Delay[J].1EEE Transaction on Acoustics, Speech and SignalProcessing, 1976, 24(4): 320-327.
发明内容
发明目的:针对现有技术中存在的问题与不足，本发明提供一种用于水声被动定位的信号到达时延差跟踪算法。在基于TDOA的水声被动定位算法处理过程中，利用相关函数全局最大值和最小值之间的幅度、时延等相对关系，辨别相关峰的真伪，从而找到真实的信号到达时延差，并对该时延差进行跟踪，输出稳定的信号到达时延差估计值。技术方案:一种用于水声被动定位的信号到达时延差跟踪算法，步骤如下:(I)变量初始化，对信号进行分帧；(2)采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函数值序列；(3)搜索互相关函数值序列的全局最大值，利用互相关函数值序列的下标与到达时延差之间的对应关系，确定当前信号帧到达两个水听器的时延差；(4)重复步骤(2)和(3)，直至得到连续5个稳定的到达时延差估计值，将最新的到达时延差估计值作为跟踪初始值；(5)使用离散傅里叶变换方法估计新一帧信号的互相关函数值序列，在跟踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值；(6)重复步骤5，直至所有信号帧处理完成。详细实施步骤为:1、步骤(I)具体为:
I)定义一个数组TraceArray含有5个元素，初值都置为O ;定义计数器变量Counter,置为O ;定义计数器变量FrameCounter,置为I ;2)设两个水听器之间的距离为d(m)，水中声速为cUXs—1)，对两个水听器接收到的水声信号的采样频率为fs(Hz)，信号长度为每个通道L点，&为极大值搜索门限，每帧数据长度为N点，N>0.1fs且N〉(d+gL) fs/c+l，帧数M为不大于L/N的最大正整数；令gff=d *fs/c,P为不大于gw的最大正整数；gP为相关函数序列全局最大值与最小值时间差的绝对值门限；3)数字滤波器H的阶数为k，下限频率为& (Hz)，上限频率为fH (Hz),且有4〈fH〈fs/4，数字滤波器采用FIR滤波器的窗函数法设计。2、步骤(2)具体为:I)设各水听器接收到的第FrameCounter巾贞信号为Sij, i表示通道序号，i=l, 2, j表示数据巾贞号，j=FrameCounter ；2)对Sij通过数字滤波器H，得到经过滤波的信号S’ u ;3)取S’ !j和S’ 2J这两个通道的信号，采用互相关函数的傅里叶变换估计法计算s’ ijPS’ 2」两个通道的互相关函数值序列R21(HI)，m为整数，且有-(N-1) ^ N-103、步骤(3)具体为:I)求取R21 (m)在下标区间[_ (N-1) +P，N-1-P]内的最大值Rmax，Rmax对应的下标为IHmax _(N~l) +P ^ HImax ^ N_1_P,对应的信号到达两个水听器的时延差为mMX.c/fs ；2)求取R21 (m)在下标区间[-(N-1)+P，N-1-P]内的最小值RMIN，Rmin对应的下标为mMiN, _ (N-1) +P ^ mMIN ^ N-1-P。
4、步骤(4)具体为:I)如果 |mMAX_mMIN| 彡 gP.fs,则 Counter 值力卩 I,并且令TraceArray (Counter) =Himax.c/fs ；2)如果 |mMX_mMIN| <gP. ；,若 IRmaxI 彡 | Rmin | ,则 Counter 值加 I,并且令 TraceArray的第 Counter 个元素 TraceArray (Counter) =HimaxC/fs ;若 I Rmax I〈 I Rmin I，则 Counter 值置为 O ;3)如果 Counter 值为 5,则进行第 5)步；如果 Counter〈5,则 FrameCounter 加 I,取第 FrameCounter 巾贞数据 Sij (i=l, 2, j=FrameCounter)；4)重复步骤(2)和(3),以及步骤(4)的第I) 3)步,直至Counter值为5 ；5)令 δ 丨=I TraceArray(2)-TraceArray (I) |, δ 2= TraceArray(3)-TraceArray (2)I，δ 3= TraceArray(4)-TraceArray (3) |，δ 4= TraceArray(5)-TraceArray (4) | ；6)若均小于 gL, Counter 置为 0,令 FrameCounter 加 I,令跟踪初始时延差为τ Q=TraceArray (5)/c,进行步骤(5)；7)若 δ 1、δ 2、δ 3、δ 4 不均小于 gL, Counter 置为 O,令 FrameCounter 加 I,返回至步骤(2)。5、步骤(5)具体为:I)设各水听器接收到的第FrameCounter巾贞信号为Sij, i表示通道序号，i=l, 2, j表示数据巾贞号，j=FrameCounter ；
2)对Sij通过数字滤波器H，得到经过滤波的信号S’ u ;3)取S’ !j和S’ 2J这两个通道的信号，采用互相关函数的傅里叶变换估计法计算s’ ijPS’ 2」两个通道的互相关函数值序列R21(HI)，m为整数，且有-(N-1) ^ m ^ N-1 ；4)令P。为不大于τ。.fs的最大正整数，Q0为不大于gL.fs/c的最大正整数，求取 R21(m)在下标区间[P0-Q0, P0+Q0]内的最大值 R’ max 及下标 m’ mx, P0-Q0 < m’ Mx < P0+Q0，输出当前巾贞时延差估计结果T=m’MX/fs,令Ttl=T ；6、步骤(6)具体为:FrameCounter加I,若FrameCounter ^ M,则返回至步骤(5);否则，结束处理。为了进一步说明本发明，现将本发明的原理介绍如下:在水声传播环境中，声源发出的信号到达两个水听器所经过的路径是不同的，传播所经过的时间延迟也是不同的。当两个水听器距离在一定的范围内，这两个水听器接收到的信号能够保持较好的相关性，通过寻求两个水听器接收到的信号的最大相似性发生的时间差，就可以得到信号到达两个水听器的时延差，这就是基于互相关函数的TDOA估计算法的基本原理。但是，声源发出的信号到达某个水听器的路径也不是单一的，有可能存在直达路径(D，即声波在传播中没有和海底或海面发生碰撞)、海面反射路径(S1，声波与海面碰撞一次后，经过反射到达水听器)、海底反射路径(B1，声波与海底碰撞一次后，经过反射到达水听器)、海面-海底反射路径(S1B1，声波先与海面碰撞一次，再与海底碰撞一次，然后到达水听器)、海底-海面反射路径(B1S1，声波先与海底碰撞一次，再与海面碰撞一次，然后到达水听器)以及多次反射路径(SnBn，声波与海底、海面发生多次碰撞后，到达水听器)等多种传播路径，通常情况下D路径的传播时延是最小的。经过这些路径传播的信号之间也具有一定的相关性，因此在采用互相关法估计TDOA时，可能这样的情况:A水听器接收到的信号中的D路径与B水听器接收信号中的B1路径相关性较好，导致互相关函数的最大值对应于B的B1路径相对于A的D路径的时延差，而真实的传播时延应该是A水听器接收到的信号中的D路径与B水听器接收信号中的D路径的时延差，这时TDOA估计就发生了错误。理论分析表明，同种性质的路径相关得到的时延差几乎相同，其估计误差可以忽略(例如A水听器接收到的信号中的B1路径与B水听器接收信号中的B1路径相关得到的时延差，与两个水听器的D路径相关得到的时延差可认为是相同的)。由于海底和海面的界面并不是平整的，经过海底反射路径和海面反射路径传播的信号幅度并不稳定，因此这种由于不同性质的路径相关得到的到达时延差通常不是稳定存在的，某些时刻比较明显，某些时刻不明显；而直达路径不和界面发生作用，经过该路径传播的信号相对比较稳定，直达路径之间相关得到的时延差是比较稳定的。鉴于水声信号传播的这种特性，可以比较连续几次的TDOA估计结果，如果这几次的估计值比较接近，则可认为是同种性质的路径相关得到的，是正确的TDOA估计值。此外，声波经过海面反射后，其相位通常会发生180度的相移，因此D路径与S1路径相关得到的值一般小于零，即对应于相关函数的最小值，根据这一点可以进行D路径与S1路径相关峰的分辨。有益效果:与现有技术相比，本发明提供的用于水声被动定位的信号到达时延差跟踪算法，能够有效地解决基于到达时延差的水声被动定位算法中的多路径干扰问题，且不需要人工指定跟踪算法的到达时延差初始值，输出的TDOA估计值稳定。
图1为本发明实施例的方法流程图；图2为水听器A接收的水声信号；图3为水听器B接收的水声信号；图4为经过带通滤波的水听器A接收信号；图5为经过带通滤波的水听器B接收信号；图6为经过带通滤波的第I帧水听器A接收信号和水听器B接收信号互相关序列，图7为经过带通滤波的第2帧水听器A接收信号和水听器B接收信号互相关序列；图8为经过跟踪处理的信号到达时延差；图9为未经过跟踪处理的信号到达时延差。
具体实施例方式下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。根据图1，介绍本实施例的用于水声被动定位的信号到达时延差跟踪算法，包括如下步骤:( I)变量初始化,对信号进行分巾贞;具体为:I)定义一个数组TraceArray含有5个元素，初值都置为O ;定义计数器变量Counter,置为O ;定义计数器变量FrameCounter,置为I ;
2)设水听器A和水听器B之间的距离d=22.7 (m)，水中声速c=1530 (mXs-1)，对两个水听器接收到的水声信号的采样频率fs=200 (kHz)，信号长度为每个通道L点，L=19064509点，极大值搜索门限gL=0.3，每帧数据长度为N点，N=40000点，帧数M=476 ；gff=2967.3, P=2967, gP=0.00025 (s)；3)数字滤波器H的阶数k=512，下限频率4=4 (kHz)，上限频率fH=12 (kHz);数字滤波器采用FIR滤波器的窗函数法设计。(2)采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函数值序列；水听器A的接收信号用Su表示，S11的波形如图2所示，水听器B的接收信号用&表示，S21的的波形如图3所示；S’ n的信号波形如图4所示，S’ 21的信号波形如图5所示；S’ n和S，21的互相关序列R21如图6所示，-39999 ^ m ^ 39999 ；(3)搜索互相关函数值序列的全局最大值，利用互相关函数值序列的下标与到达时延差之间的对应关系，确定当前信号帧到达两个水听器的时延差；Rmx=475.1，mMX=1842,Rmin=-484.6, mMIN=275 ；(4)重复步骤(2)和(3)，直至得到连续5个稳定的到达时延差估计值，将最新的到达时延差估计值作为跟踪初 始值；mMX-mMIN| =1570>gP.fs=50,则 Counter 值加 I, TraceArray (I) =mMX.c/fs=14.0913 ；当 Counter 值为 5 时，TraceArray (2) =-14.0925, TraceArray (3) =-14.0886,TraceArray(4)=-14.0896, TraceArray(5)=-14.0949, δ ^0.0012, δ2=0.0039, δ 3=0.001,δ 4=0.0053 ； δ ^ δ 2> δ 3> δ 4 均小于 gL,贝丨J Counter 置为 0, FrameCounter 加 I，令跟踪初始时延差为Ttl=-0.009212。此处到达时延差为负值，表示信号到达B水听器的时间早于到达A水听器的时间；(5)使用离散傅里叶变换方法估计新一帧信号的互相关函数值序列，在跟踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值；S，12 和 S，22 的互相关函数值序列 R21 如图 7 所示，P0=1842, Q0=39, R，Mx=524.4,m，Mx= 1842, R，MIN=-460.0, m，MIN=1826, τ =m’ mx/fs=-0.009212,令 τ 0= τ ；(6)重复步骤5，直至所有信号帧处理完成；FrameCounter 加 I,若 FrameCounter ^ 476,则返回至步骤(5);否则，结束处理。最后得到的信号到达时延差序列如图8所示，此处为了显示跟踪处理的效果，将未经跟踪处理的信号到达时延差估计结果也一并画出，如图9所示。
权利要求
1.一种用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤如下: (1)变量初始化，对信号进行分帧； (2)采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函数值序列； (3)搜索互相关函数值序列的全局最大值，利用互相关函数值序列的下标与到达时延差之间的对应关系，确定当前信号帧到达两个水听器的时延差； (4)重复步骤(2)和(3)，直至得到连续5个稳定的到达时延差估计值，将最新的到达时延差估计值作为跟踪初始值； (5)使用离散傅里叶变换方法估计新一帧信号的互相关函数值序列，在跟踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值； (6)重复步骤5，直至所有信号帧处理完成。
2.如权利要求1所述的用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤(I)具体为: 1)定义一个数组TraceArray含有5个元素，初值都置为O;定义计数器变量Counter，置为O ;定义计数器变量FrameCounter，置为I ; 2)设两个水听器之间的距离为d(m)，水中声速为c (mX s—1)，对两个水听器接收到的水声信号的采样频率为fs(Hz)，信号长度为每个通道L点，&为极大值搜索门限，每帧数据长度为N点，N>0.14且吣((^^4/(3+1，帧数1为不大于1^的最大正整数；令8 =(1.^/c,P为不大于gw的最大正整数；gP为相关函数序列全局最大值与最小值时间差的绝对值门限；` 3)数字滤波器H的阶数为k，下限频率为fjHz)，上限频率为fH(Hz)，且有4〈fH〈fs/4，数字滤波器采用FIR滤波器的窗函数法设计。
3.如权利要求2所述的用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤(2)具体为: 1)设各水听器接收到的第FrameCounter巾贞信号为Sij,i表示通道序号，i=l, 2, j表示数据巾贞号，j=FrameCounter ； 2)对Sij通过数字滤波器H，得到经过滤波的信号S’u ; 3)取S’!j和S’ 2J这两个通道的信号，采用互相关函数的傅里叶变换估计法计算S’ !j和S’ 2J两个通道的互相关函数值序列R21 (m), m为整数,且有-(N-1)彡m彡N-1。
4.如权利要求3所述的用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤(3)具体为: 1)求取R21(Hi)在下标区间[-(N-l)+P，N-1-P]内的最大值Rmx，Rmax对应的下标为IHmax _(N~l) +P ^ HImax ^ N_1_P,对应的信号到达两个水听器的时延差为mMX.c/fs ； 2)求取R21(Hi)在下标区间[-(N-l)+P，N-1-P]内的最小值Rmin，Rmin对应的下标为mMiN, _ (N-1) +P ^ mMIN ^ N-1-P。
5.如权利要求4所述的用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤(4)具体为:I)如果 I mMX-mMIN I 彡 gP.fs,则 Counter 值加 I,并且令 TraceArray (Counter) =mMX.c/fs；.2 )如果 I HImax-HImin I <gP.fs,若 I Rmax I ≥ I Rmin I，则 Counter 值加 I,并且令 TraceArray 的第 Counter 个元素 TraceArray (Counter) =mMxc/fs ;若 I Rmax I〈 I Rmin I，则 Counter 值置为 O ; .3)如果Counter值为5,则进行第5)步；如果Counter〈5,则FrameCounter加I,取第FrameCounter 巾贞数据 Sij (i=l, 2, j=FrameCounter)； .4)重复步骤(2)和(3),以及步骤(4)的第I) 3)步,直至Counter值为5；.5)令S1=ITraceArrayO-TraceArray(I) , δ 2= | TraceArray (3) -TraceArray (2) |，δ.3=I TraceArray (4)-TraceArray (3) |，δ 4= TraceArray(5)-TraceArray (4) | ； .6)若δρ δ 2、δ 3、δ 4均小于gL, Counter置为0,令FrameCounter加I,令跟踪初始时延差为 τ Q=TraceArray (5)/c,进行步骤(5)； .7)若δ^ δ δ δ 4不均小于gL, Counter置为O,令FrameCounter加I,返回至步骤(2)。
6.如权利要求5所述的用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤(5)具体为: .1)设各水听器接收到的第FrameCounter巾贞信号为Sij,i表示通道序号，i=l, 2, j表示数据巾贞号，j=FrameCounter ； .2)对Sij通过数字滤波器H，得到经过滤波的信号S’u ; .3)取S’!j和S’ 2J这两个通道的信号，采用互相关函数的傅里叶变换估计法计算S’ !j和S’ 2J两个通道的互相关函数值序列R21 (m), m为整数,且有-(N-1) ≤ m ≤ N-1 ； .4)令Ptl为不大于τ^ 的最大正整数,Qtl为不大于*fs/c的最大正整数,求取R21 (m)在下标区间[P0-Q0，VQ0]内的最大值R， Αχ及下标m’祖，P0-Q0≤m’ Mx≤Pq+Qq，输出当前中贞时延差估计结果τ =m’ MAX/fs,令τ f τ。
7.如权利要求6所述的用于水声被动定位的信号到达时延差跟踪算法，其特征在于步骤(6)具体为: FrameCounter加I,若FrameCounter ≤ Μ,则返回至步骤(5);否则，结束处理。
全文摘要
本发明公开了一种用于水声被动定位的信号到达时延差跟踪算法，能够有效地解决基于到达时延差的水声被动定位算法中的多路径干扰问题，且不需要人工指定跟踪算法的到达时延差初始值。包括如下步骤(1)变量初始化，对信号进行分帧；(2)采用离散傅里叶变换方法估计两个水听器接收到的信号帧的互相关函数值序列；(3)确定当前信号帧到达两个水听器的时延差；(4)将最新的到达时延差估计值作为跟踪初始值；(5)估计新一帧信号的互相关函数值序列，在跟踪初始值对应的相关函数值序列下标的附近搜索互相关函数值序列的局部最大值并计算到达时延差，作为新一帧信号到达时延差的估计值，并更新跟踪初始值，直至所有信号帧处理完成。
文档编号G01S5/18GK103176166SQ20131008064
公开日2013年6月26日 申请日期2013年3月14日 优先权日2013年3月14日
发明者安良 申请人:东南大学