强跟踪滤波滤波器

强跟踪滤波滤波器（共8篇）

强跟踪滤波滤波器 篇1

对活动目标进行跟踪是计算机视觉领域一个非常重要的研究方向,目标跟踪的实质就是对目标状态的估计、预测和滤波,广泛应用于视频监控系统、飞行控制系统及来袭目标检测等领域。大多数情况下,对匀速运动目标要求具有良好的跟踪能力,同时也要求对复杂环境和机动目标也要有较强的跟踪能力。

早期的视频跟踪系统中,运用经典卡尔曼滤波器(Kalman Filter,KF)即无偏最小方差估计,可以得到较好的跟踪效果。然而,KF关于模型和噪声统计特性不确定性的鲁棒性比较差,当视频图像序列中的运动目标发生突变时,KF进行跟踪可能会出现运动目标丢失现象,导致滤波误差增大、跟踪性能下降。为此,笔者提出将强跟踪滤波器应用于匀速运动目标和机动性目标的跟踪,并对其进行仿真。

1 KFKF是一种状态预测估计算法,假设一个运动系统的状态方程为:

$\{\begin{array}{l}X(k+1)=F(k)X(k)+V(k)\\ z(k+1)={\rm H}(k)X(k+1)+W(k+1)\end{array}$

(1)

采用KF跟踪其运动轨迹的方法可推导为:

$\begin{array}{l}\widehat{X}(k|k+1)=F(k)\widehat{X}(k|k)\\ {\rm P}(k+1|k)=F(k){\rm P}(k|k)F{}^{\text{Τ}}(k)+Q(k)\\ {\rm K}(k+1)={\rm P}(k+1|k){\rm H}{}^{\text{Τ}}(k)[{\rm H}(k){\rm P}(k+1)\cdot \\ {\rm H}{}^{\text{Τ}}(k)+R(k+1)]{}^{-1}\\ y(k+1)=z(k+1)-{\rm H}(k)\widehat{X}(k+1|k)\\ \widehat{X}(k+1|k+1)=\widehat{X}(k+1|k)+{\rm K}(k+1)\gamma (k+1)(2)\\ {\rm P}(k+1|k+1)=[{\rm I}-{\rm K}(k+1){\rm H}]{\rm P}(k+1|k)\end{array}$

$\{\begin{array}{l}R(k)=E\{W(k)W{}^{\text{Τ}}(k)\}\\ Q(k)=E\{V(k)V{}^{\text{Τ}}(k)\}\end{array}$

式(1)所示的数学模型具有很大的不确定性,其原因主要是:模型过于简化,系统噪声的统计特性不是很准确以及对现实系统初始状态的统计特性的建模较粗糙以及现实系统的参数随时间发生变动等。KF对于模型不确定性的鲁棒性很差,导致状态预测估计不准,甚至出现错误发散。因此,KF只适用于符合高斯、线性化的情况。因此又推出了扩展KF。

扩展KF的原理是,当系统达到平稳状态时,系统会丧失对突变状态的跟踪能力,这制约了扩展KF(也包括KF)在视频目标中的跟踪应用。究其原因,是由于系统达到平稳状态时,扩展KF的增益矩阵K(k+1)将趋于极小,如果此时系统状态发生突变,预报残差γ(k+1)将随突变增大,但其增益矩阵K(k+1)仍将保持为极小值。因此,从式(2)推导出的扩展KF将丧失对突变状态的跟踪能力。为此,可以认为扩展KF是一种非闭环滤波器,原因是这种类型的滤波器的K(k+1)不会随着滤波效果自适应地进行改变,以确保始终保持系统状态的准确跟踪能力。这一现象在线性定常随机系统中更为突出。此时,只需用普通KF(或扩展KF)进行状态预测估计。而KF的增益矩阵可以根据线性系统的参数阵(A,B,C)离线计算,然后存储在计算机内存中被在线应用,这时如果系统状态发生突变,滤波器的增益矩阵不会随之改变(事先存好的静态变量),因此KF也就丧失了对突然变化状态的跟综能力,故扩展KF也是一种非闭环滤波器。

2 强跟踪滤波器

如果一个滤波器与通常的滤波器相比满足:很好的关于模型不确定性的鲁棒性;极强的关于突变状态的跟踪能力,即使在系统达平稳状态时仍保持对缓慢变化状态与突变状态的跟踪能力;合适的计算复杂性。那么就可以确定该滤波器为强跟踪滤波器(Strong Tracking Filter,STF)。前两个条件就是为了解决扩展类KF的两大问题提出来的,最后一个条件则是为了保证STF的实际应用。

由式(1)确定的一类系统的STF方程应具有下述的一般结构:

$\begin{array}{l}\widehat{X}(k+1|k+1)=\widehat{X}(k+1|k)+{\rm K}(k+1)\gamma (k+1)\\ \widehat{X}(k|k+1)=f[k,u(k),\widehat{X}(k|k)]\\ \gamma (k+1)=y(k+1)-h[k+1,\widehat{X}(k+1|k)]\end{array}$

现在需要解决的难题是,要实时在线地确定时变增益矩阵K(k+1)以使此滤波器满足STF的所有条件。为此,提出正交性原理,即为了使滤波器(2)为STF的一个充分条件是在线选择一个适当的实时变化的增益矩阵K(k+1),使得:

$\begin{array}{l}E\{[X(k+1)-\widehat{X}(k+1|k+1)][\widehat{X}(k+1)-\widehat{X}(k+1|k+1){}^{\text{Τ}}]\}=\min (3)\\ E\{\gamma (k+1+j)\gamma {}^{\text{Τ}}(k+1)\}=0,k=0,1,2,\cdots ；\end{array}$

j=1,2,… (4)

条件(4)要求在不同时刻的残差序列处处相互正交。条件(3)本来就是原扩展类KF的性能指标。

现说明此正交性原理的一个例子:已经证明,当参考模型与实际系统完全匹配时,KF输出的残差序列是一系列的非自相关的高斯白噪声序列。所以,式(3)是满足的。而式(4)本来就是KF的性能指标,因此也是满足的。

由于所建模型不确定性的影响,使得滤波器的状态估计值在偏离系统状态时,一定会在输出的残差序列的均值和幅值上表现出来。在这时,如果在线实时地调整增益矩阵,强迫式(4)仍然成立,造成残差序列仍然保持相互正交,就可以强迫STF保持对实际系统状态的跟踪。

正交性原理拥有很强的物理意义,它说明当存在模型的不确定性时,应在线实时调整增益矩阵K(k+1),迫使输出的残差序列始终保持具有类似高斯白噪声的性质。这也说明已经将输出的残差序列中的所有有效信息都提取出来了。

在不存在所建模型的不确定性时,STF按照正常情况运行,式(4)已经满足,起不到调节作用。此时STF就演化为一般的基于性能指标(3)的扩展类KF。

这类正交性原理的核心是式(4)和当用其他的性能指标取代式(3)后,就能够得到另外一些变了形的类似的正交性原理。所以,当在原先的滤波器上额外加上条件(4)后,就可以改变原先的滤波器,使得其满足STF的所有条件。

对于非线性时变系统,在实际中应用该正交性原理时,式(3)、(4)非常难精确地满足。这时要使其近似满足即可,以便减少计算量,使强跟踪滤波器能够实时在线进行视频跟踪。

3 仿真

为了验证上述算法的有效性,进行仿真分析。试验用的处理器为主频3.00GHz Intel(R)E8400 CPU,内存为2GHz的PC机;所有算法由Matlab 7.1编程实现。

由运动模型产生的数据进行仿真分析,能够更好地看出跟踪效果并进行分析。

图1、2是分别使用KF和STF在不同情况下的仿真结果。

由图1b、c可以看出,在目标运动变化较缓慢时,KF和STF对目标的跟踪效果没有太大的差别。由图2b、c可以看出,在观测数据受到严重污染、干扰和目标大机动运动时,STF的跟踪效果明显好于KF的跟踪效果,在实际的目标跟踪中具有很强的优势。在目标发生大的机动变化时,STF仍然具有较好的目标跟踪能力,这就是STF具有而KF或扩展类KF所不具有的优点。其中STF的卡尔曼增益矩阵在运算时,当系统出现大的机动变化时取无穷大,以便对运动系统的大变化情况进行跟踪,之后KF的卡尔曼增益又恢复到极小值。

4 结束语

仿真结果表明,强跟踪滤波器在目标跟踪应用中比卡尔曼滤波器具有较强的鲁棒性和对噪声信号的适应性,能够更好地解决卡尔曼滤波器和扩展类卡尔曼滤波器所遇到的问题。

摘要：针对早期卡尔曼滤波在视频跟踪中的应用情况及其存在的缺陷,将强跟踪滤波器用于匀速运动目标和机动性目标的跟踪。仿真结果表明:在一定条件下,强跟踪滤波器具有较强的抗干扰能力和跟踪性能。

关键词：强跟踪滤波滤波器,匀速和机动性目标,视频跟踪,仿真

强跟踪滤波滤波器 篇2

粒子滤波器的优化及在纯方位跟踪中的应用

为了克服粒子滤波在应用中由于建议分布函数选择不合理导致粒子数减少进而使其状态量失去多样性的问题,提出了改进建议分布的混合粒子滤波方法(UKPF).首先,在粒子滤波的基础上融合进扩展卡尔曼滤波及无迹卡尔曼滤波;然后,融合后的新算法在计算建议概率密度分布时,粒子的.产生充分考虑当前时刻的量测,利用所提混合算法来加入最新的观测量并产生粒子滤波的建议分布;最后,使得粒子的分布更加接近状态的后验概率密度,最大化地实现其滤波性能.仿真结果表明:对于纯方位跟踪问题,所提算法不仅解决了EKF的线性化损失问题及UKF在解决一般非高斯问题中建模的困难,而且与PF等粒子滤波器相比,具有更高的跟踪精度.

作 者：许立太 作者单位：兰州石化职业技术学院,甘肃,兰州,730060刊 名：探测与控制学报 ISTIC PKU英文刊名：JOURNAL OF DETECTION & CONTROL年，卷(期)：32(2)分类号：V241关键词：粒子滤波 建议分布 纯方位跟踪 机动目标

强跟踪滤波滤波器 篇3

在编队卫星相对导航的研究中扩展卡尔曼滤波(extended Kalman filter,EKF)因算法简单而得到广泛的应用[5,6],但因其对非线性函数泰勒展开式进行一阶线性化截断而降低了滤波估计精度,并且在系统模型存在不确定性时易于发散,鲁棒性不强。进而,文献[7,8]分别研究了编队卫星相对导航的鲁棒EKF算法和自适应EKF算法。为了克服EKF线性化截断的缺陷,一些研究中采用Sigma点卡尔曼滤波(sigma point Kalman filter,SPKF)进行相对导航估计,包括基于UT(unscented transform)变换的无迹卡尔曼滤波(unscented Kalman filter,UKF)[9,10]和基于Stirling多项式插值的中心差分卡尔曼滤波(central difference Kalman filter,CDKF)[11],它用Sigma点的分布近似表示非线性函数的分布,有效地提高了编队卫星相对导航估计精度,但是SPKF算法仍然要求精确已知噪声的先验统计特性。然而,在实际应用中,因系统受不确定因素的影响,噪声的先验统计特性极易发生变化。为提高相对导航滤波算法应对系统模型不确定性及噪声统计特性不准确的能力,鲁棒H∞滤波算法[12]、容积卡尔曼滤波(cubature Kalman filter,CKF)[13]算法和自适应容积卡尔曼滤波(adaptive cubature Kalman filter,ACKF)[14]算法分别被应用于编队卫星相对导航的研究中。作为提高算法鲁棒性的一种重要方法,强跟踪滤波(strong tracking filer,STF)[15]算法通过在状态预测协方差阵中引入渐消因子来在线实时调整增益矩阵,强迫输出残差序列保持相互正交,从而实现在系统模型不确定时对系统状态的跟踪。但是STF要求非线性函数连续可微以计算其Jacobi矩阵并且对非线性状态后验分布近似精度只能达到一阶。因此,将STF算法与UKF算法相结合,可以在保持滤波算法鲁棒性的同时有效提高滤波估计精度,减少计算量并拓展滤波算法应用范畴。结合以上分析,本文提出一种简化强跟踪UKF(SS-TUKF)滤波算法。

1 自主相对导航数学模型

1.1 编队卫星动力学模型

主星Hill坐标系定义为:坐标系原点O位于主星的质心,Ox轴由质心指向从星,Oy轴在轨道面内垂直于Ox轴,指向主星运动方向,Oz轴与Ox轴、Oy轴构成右手坐标系,如图1所示。

在主星Hill坐标系中,定义从星相对轨道状态,则相对动力学方程为[16]:

式中,θ为主星真近点角;r为主星地心距;μ为地球引力常数。考虑主星运行于近圆轨道,则,将式(1)右端在原点处泰勒展开,保留一阶项,整理得到描述卫星相对运动的CW方程如下:

将CW方程写成离散线性形式为:

式中,wk-1为未建模的摄动加速度,可以等效成零均值的高斯白噪声;G为如下状态转移矩阵:

1.2 自主量测模型

使用激光测距仪和光学相机等设备可以直接测量从星与主星之间的视线距离ρ,方位角φ和俯仰角θ,如图2所示,测量值与相对轨道状态x有如下非线性关系:

写成量测方程的形式如下:

式(8)中,v为量测噪声。

2 相对导航滤波算法

传统EKF算法对模型不确定性的鲁棒性较差,而STF对模型参数失配具有强鲁棒性并且对系统量测噪声时变统计特性具有低敏感性。实现STF的关键是在线实时调整增益矩阵Kk使得残差序列相互正交,即满足如下正交性原理:

式中,条件1)是基本的性能指标,条件2)要求输出残差序列时刻保持正交。

2.1 标准强跟踪滤波算法

STF算法可以在EKF基础上通过在状态误差协方差阵中引入时变渐消因子来实现。已知系统状态方程和量测方程如下:

标准强跟踪滤波算法过程如下。

步骤1滤波器的初始化。

循环k=1,2,…,完成以下步骤。

步骤2计算如下雅克比矩阵。

步骤3引入渐消因子进行时间更新。

在估计误差协方差矩阵中引入次优渐消因子,则状态估计和估计误差协方差时间更新为:

式中,渐消因子按如下方法计算。

首先,估算残差协方差矩阵Vk=E(εkεkΤ)为:

式(19)中,0<ρ≤1为遗忘因子,通常取ρ=0.95。

然后,定义矩阵Mk,Nk如下:

式中,β≥1为弱化因子,其意义是避免渐消因子引起过调节作用,使状态估计更加平滑,可以由以下准则借助计算机仿真确定:

最后,计算渐消因子λk的次优解为:

步骤4状态量测更新。

按如下滤波方程进行状态估计的量测更新和估计误差协方差的更新:

2.2 简化强跟踪UKF滤波算法

STF算法对非线性状态后验分布近似精度只有一阶且并需计算Jacobi矩阵。理论上已经证明,UT变换无需计算Jacobi矩阵,对后验均值和协方差的近似精度可以达到三阶,且适用于不连续可微的非线性系统。因此,采用UT变换代替Jacobi矩阵的求解可以得到基于UT变换的强跟踪UKF滤波算法。然而,在强跟踪UKF算法中,无论状态方程是否为线性均对其进行UT变换,而本文研究的编队卫星相对动力学方程为线性方程,对线性系统而言,卡尔曼滤波可以得到理论上的最优估计,此时若进行无迹变换反而会降低估计精度,因此可以用标准卡尔曼滤波中的时间更新代替UKF中相应的无迹变换过程,这样可以有效地降低运算量并提高估计精度,增强算法的实时性。

为了推导简化强跟踪UKF滤波算法,首先需要进行STF的等价表示,Pk,Pz,Pxz可以表示如下:

式中,上标(l)表示没有引入渐消因子时的相应的量。前面已经假设Qk-1为正定矩阵,故Pk-(l)的逆矩阵一定存在,由式(29)得:

代入STF的Nk,Mk表达式得:

到此,便得到了Mk、Nk的等价表达式,次优渐消因子λk的计算方法同式(23)。

已知n维离散时间非线性系统如下:

简化强跟踪UKF滤波算法如下。

步骤1滤波器的初始化。

步骤2无渐消因子k时刻时间更新。

用标准卡尔曼滤波计算状态预测值与预测误差协方差矩阵如下:

步骤3计算渐消因子λk及带渐消因子的Pk-。

首先根据xk的均值和协方差当前最好估计是和Pk-,选择sigma点如下:

将sigma点转换为^zk(i)(l)如下:

合并得到k时刻量测预测和协方差预测为:

其权值

估计之间的协方差为:

到此可以计算λk,为使滤波器具有对系统模型不确定性的鲁棒性,在状态预测协方差阵Pk-(l)中引入渐消因子λk,进而可以计算Pk-为:

步骤4计算协方差矩阵Pz和Pxz。

利用带渐消因子的Pk-重新计算协方差矩阵Pz和Pxz如下所示。

步骤5状态量测更新。按如下滤波方程进行状态的量测更新:

3 仿真试验结果分析

编队卫星中主星和从星的轨道根数如表1所示,此时系统模型中的n=0.001 038,用STK中HPOP预报器生成高精度相对轨道并将其作为真实值,建模时考虑地球非球形摄动、大气阻力摄动和三体引力摄动的影响。假设测距精度为1 cm,测角精度为0.01°,采样周期为1 s。相对导航初始位置和速度的真实值分别为[-81.301 116.672 140.641]Τ和[0.060357 0.168568-0.104951]Τ,假设主星不能准确获得从星的初始相对位置和速度,即假设状态估计初值

,仿真时间为3 000 s,同时采用EKF和SSTUKF两种方法进行对比分析。

图3、图4给出了整个仿真周期内相对位置和相对速度的估计误差曲线,从图中可以看出EKF和SSTUKF两种方法均可使相对状态估计误差收敛到零附近,从而实现相对导航。将仿真前20 s和后80s放大,如图5、图6所示。从图中可以看出,EKF的收敛速度与SSTUKF基本相当,均在10 s左右收敛,但滤波估计精度低于SSTUKF。统计2 500~3 000 s的估计误差均方差如表2所示,从表中可以看出,相比于EKF,SSTUKF的估计精度大约提高了25%。

为了验证SSTUKF算法在模型参数失配时的鲁棒性以及对时变噪声统计特性的低敏感性,将系统模型中的n值改为0.02,此时的系统模型参数与真实的系统模型参数已有较大的偏差,同时将噪声标准差乘以0.2sin(0.01t)以使其具有时变性质,将此条件称为情况二(相应的,前述参数正常时称为情况一),其他参数不变,仿真结果如图7和图8所示。统计2 500 s~3 000 s的估计误差均方差如表3所示,结合图表可以看出,在模型参数偏差较大及噪声时变条件下,情况二的滤波估计精度仅比情况一略低,从而验证了SSTUKF相应的鲁棒性与低敏感性。

4 结论

雷达跟踪滤波技术研究 篇4

雷达数据处理技术最早可追溯至19世纪初高斯提出的最小二乘法，其运用最小二乘法开创了用数学方法处理观测和实验数据的科学领域。之后对这种方法进行了不断的修改和完善，成为了现代滤波理论的基础。在现代雷达系统中，信号处理的检测解决的是有无目标存在的问题，而数据处理则是在目标已经存在的前提下，给出其参数估计。

作为雷达信号处理的后续处理过程，雷达数据处理将数据关联技术与现代滤波理论有机结合用于消除由杂波剩余等造成的虚假检测，改善雷达信号处理性能;根据雷达的测量值利用参数估值理论从大量雷达回波点迹中提取捕获目标点迹，进行点迹与点迹关联，进行目标航迹的建立及管理，确定目标的位置、速度、机动情况等参数;进行点迹与航迹关联、跟踪滤波、航迹维持，形成综合目标态势：进行目标航迹预测外推，调度系统资源对多目标进行跟踪监视，形成连续的目标航迹数据，向火控和武器系统输出高数据率目标指示信息。

2 跟踪滤波算法

跟踪滤波是雷达数据处理技术的重要组成部分，其主要功能是根据目标的实际测量信息实时估计目标当前的位置、速度等参数，并外推出下一次天线扫描周期目标出现的位置信息。该外推信息在跟踪雷达中用于检验下一拍测量信息的合理性;在搜索雷达中用于航迹的相关处理。常见的跟踪滤波器有α-β滤波器、Kalman滤波器等，根据不同的计算资源、需处理的目标数量、目标的不同动态特性、雷达精度要求等条件可选择不同的跟踪滤波器。

与Kalman滤波器相比，α-β滤波器的优点是算法较简单、计算工作量较小、容易实现。在本文应用中选用α-β滤波器，其对跟踪目标进行最佳线性无偏估计，并按航迹质量和预测误差改变滤波器的参数，实现自适应滤波。在混合坐标系下建立数学模型如下：

滤波方程：

其中:k为雷达的观测次数

滤波的外推方程:

的获取采用工程上常用的经验公式:

其中，r为信噪比，C为大于零的常数。

3 典型航迹仿真

根据跟踪滤波算法建立的数学模型，选取几种典型的雷达模拟航迹在Matlab工具中进行仿真，进一步分析该滤波算法的精度。

选取的典型模拟航迹有：水平匀速直线运动目标、水平蛇形机动目标、高度方向上跃起俯冲运动目标。将得到的目标模拟航迹数据，通过直球坐标转换后在距离R、方位角A、俯仰角E上的数据中加入测量噪声作为滤波器的输入，测量噪声的标准差分别为距离40米、方位角0.04弧度、俯仰角0.01弧度。

3.1 水平匀速直线运动目标仿真分析

选取水平匀速直线运动目标的参数为：起始位置为(x0, y0, z0)=(60km, 60km, 8km), X、Y方向的速度均为500m/s，加速度为0，航路捷径为0，采样周期为2s。

水平匀速直线运动目标的俯仰角滤波曲线如图1、距离和方位角滤波曲线如图2、距离方位角俯仰角的误差标准差曲线图3。其中图3中的蓝色曲线为测量值与真实值之差，红色曲线为滤波值与真实值之差，从中可以看出，滤波效果明显，距离误差峰值、方位角误差峰值、俯仰角误差峰值基本减小了一半以上。

3.2 水平蛇形机动目标仿真分析

选取水平蛇形机动目标的参数为：起始位置为(x0, y0, z0)=(30km, 30km, 2km)，转弯时角速度为0.1rad/s^2(相当于机动13g)，半径为3000m，采样周期为2s。

水平蛇形机动是水平机动的一种，一般发生在近程和低空，而且比水平机动的机动性更大。水平蛇形机动目标的俯仰角滤波曲线如图4、距离和方位角滤波曲线如图5、距离方位角俯仰角的误差标准差曲线图6。从图6中看出，由于水平蛇形机动目标的机动更大，距离滤波误差较大，对方位和俯仰的滤波效果明显。

3.3 高度方向上跃起俯冲运动目标仿真分析

选取高度方向上跃起俯冲运动目标的参数为：初始位置为(x0, y0, z0)=(8km, 8km, 1km)，利用峰峰高度2000，跨度100的正弦曲线模拟目标运动轨迹，最大机动约为10.5g。

高度方向上跃起俯冲运动目标的俯仰角滤波曲线如图7、距离和方位角滤波曲线如图8、距离方位角俯仰角的误差标准差曲线图9。由于目标作高度方向上的跃起俯冲运动，水平不机动，故在图7～9中距离和方位的滤波效果接近于水平匀速运动的目标。目标高度上机动和目标水平面上机动时，滤波算法对俯仰角滤波的效果都可以达到令人比较满意的效果。

3.4 典型目标的误差分析

运用蒙特卡洛方法，对上述典型目标各进行了50次仿真，将得到的误差标准差平均值整理得到表1。从中可以看出进入稳定跟踪后，滤波得到的距离值的误差标准差基本在20左右，方位角的误差标准差基本在0.015左右，俯仰角的误差标准差基本在0.0035左右。仿真结果证明，采用上述跟踪滤波算法的精度在可接受范围之内，满足数据处理的跟踪要求。

4 结语

本文介绍了雷达数据处理技术中跟踪滤波算法的基本原理，在跟踪滤波算法中选用α-β滤波器建立其在混合坐标系下的数学模型，用于实现自适应滤波;选取了典型的雷达目标模拟航迹在Matlab工具中进行仿真分析，验证了该跟踪滤波算法的精度在可接受范围之内，满足雷达数据处理的跟踪要求。

摘要：本文介绍了雷达数据处理技术中跟踪滤波算法的基本原理, 以在混合坐标系下的α-β滤波模型为研究对象, 选取了典型的三维目标进行数据建模, 验证了该跟踪滤波算法的精度满足雷达数据处理的跟踪要求。

关键词：雷达数据处理,跟踪滤波,仿真

参考文献

[1]李媛媛, 薛媛等.基于Matlab的雷达系统仿真.电子元器件应用, 2009.9.

[2]毛滔, 李盾等.相控阵雷达数据处理仿真研究.航天电子对抗, 2004 (4) .

[3]陈强超, 单雷达航迹滤波与卡尔曼滤波算法.信息化研究, 2010.5.

[4]陈明燕, 张伟等.相控阵雷达数据处理系统的仿真.电视技术, 2008.4.

强跟踪滤波滤波器 篇5

在检测过程中, 可以利用目标航迹的特征来提高系统的检测性能, 即所谓的检测前跟踪 ( TrackBefore Detect, TBD) 方法。常见的TBD方法包括Hough变换[1]、动态规划[2]及最大似然估计[3]等, 这些方法基于多帧数据沿目标的运动方向对信号进行积累, 所需的计算量较大, 且不能以迭代的方式进行处理, 很难处理复杂的目标运动模型。贝叶斯检测前跟踪算法直接对观测数据进行处理, 计算目标在各分辨单元存在的后验概率密度函数, 在这些后验概率密度函数的基础上, 对目标进行检测和跟踪, 能够以迭代的方式进行计算, 且能够处理任意的目标运动模型[4,5]。然而, 在雷达等系统的TBD应用中, 大多数文献考虑目标观测模型的方式同先检测后跟踪系统中目标的观测模型一样, 仅针对单个的分辨单元对目标的观测建立相应的概率模型[6,7], 而没有根据观测信号产生的机理建立目标的观测模型。文献[8]在雷达等系统的TBD应用中考虑了扫描波束对观测数据的影响, 即有限宽度的扫描波束能够在多个波束范围内产生目标回波, 接着应用Viterbi方法对目标进行检测。在实际的雷达等系统中, 观测数据同波束方向图及其扫描方式和匹配滤波器的时频域特性有关[9,10,11]。考虑匹配滤波器的时频域特性, 目标会在多个距离—多普勒单元产生回波, 其回波大小与信号带宽及匹配滤波器的特性有关。即一个点目标可能会分布在多个距离—多普勒分辨单元内, 在这些分辨单元内的回波也都是目标回波能量的一部分, 利用这些回波能量及对应分辨单元位置的关系可以提高系统的检测和跟踪性能, 也使得目标的观测模型与实际的观测数据更加吻合, 从而避免由于模型失配所带来的性能损失。

在贝叶斯检测前跟踪算法中, 考虑雷达等系统中匹配滤波器时频域特性对观测数据产生的影响, 并建立相应的模型, 以使目标的观测模型更符合实际观测数据, 避免由于模型失配所带来的性能损失;利用目标在多个距离—多普勒分辨单元扩散的能量及对应位置的相互关系能够提高系统的检测和跟踪性能。下面首先对贝叶斯检测前跟踪算法进行描述, 之后建立目标的状态转移方程及目标的观测模型, 其中对目标的观测模型进行详细叙述, 最后进行仿真验证。

1 贝叶斯检测前跟踪算法

利用贝叶斯理论对目标进行检测和跟踪, 其主要过程是基于观测到的每一帧数据、目标的运动模型、噪声和目标的统计特性计算每一分辨单元存在目标的后验概率密度函数, 在这些后验概率密度函数的基础上对目标检测和跟踪。贝叶斯检测前跟踪算法能够在帧数据间以概率的形式对信号进行积累, 从而可以对低信噪比的目标进行检测和跟踪[4]。

1. 1 系统模型

首先定义目标的状态矢量xk, 下标k表示采样时刻, 目标的状态转移方程可以表示为:

式中, fk - 1 ( ·) 表示相邻时刻目标状态转移的函数关系; vk表示k时刻目标动态模型的噪声, 用来衡量在两相邻时刻目标状态转移的不确定性。

在雷达等系统中, 目标回波大小是未知的, 在检测和跟踪过程中需要将目标回波的幅度Ik作为目标状态的一个分量。为了实现跟踪和检测过程的一体化处理, 需要在状态中引入一个分量Ek用来表示目标是否存在的情况 ( Ek= H1表示目标存在的情况; Ek= H0表示目标不存在的情况; H1表示目标存在的假设; H0表示目标不存在的假设) 。由此, 目标的状态矢量可以表示为xk= [mTk, Ik, Ek]T, 其中mk表示k时刻目标状态的运动分量 ( 在直角坐标系中, 当不考虑加速度时, mk=[xk, yk, zk, x·k, y·k, z·k]T, 其中mk中前3项表示目标在空间中的位置, 后3项表示目标对应的速度分量; 考虑加速度时, mk中还需加入对应的加速度分量) , 上标T表示向量或矩阵的转置。

在雷达等系统中, 针对贝叶斯检测前跟踪算法, 由于直接对原始观测数据进行处理, 处理的观测量为整个观测空域的一个多维观测数据矩阵, 下面对距离维和多普勒维进行考虑, 观测数据zk= { zrk, d} ( r, d分别为距离分辨单元和多普勒分辨单元的索引) 。由此, 对于某一分辨单元, 目标的观测模型为:

式中, wrk, d表示分辨单元的观测噪声; hk (·) 表示目标状态同观测值之间的关系。

1. 2 预测和滤波

为了以一种迭代的方式计算目标状态的后验概率密度函数, 可以分预测和滤波来进行。

预测是根据目标现在的状态, 以及目标状态的转移方程预测目标下一时刻的状态, 其计算过程为[4]:

目标状态的转移模型设定为一阶马尔可夫模型, 即目标在下一时刻的状态只与目标当前时刻的状态有关。设定系统噪声和观测噪声相互独立, 式 ( 3) 可以变为:

式 ( 3) 和式 ( 4) 中, p (xk- 1Zk- 1) 为上一时刻目标状态的后验概率密度函数; Zk - 1= ( zk - 1, …, z1) 表示直到k -1时刻所有的观测数据; p ( xk| xk - 1) 为目标状态转移概率密度函数, 根据式 ( 1) 进行求解。

滤波主要是根据系统当前时刻的观测值, 对预测值进行修正。其计算过程为[8]:

从式 ( 2) 可以看出, 目标当前时刻的观测只与目标当前时刻的状态有关。因此, 式 ( 5) 可以写为:

1. 3 检测和跟踪

在得到目标状态的后验概率密度函数之后, 就可以利用它来对目标进行检测和跟踪。在贝叶斯理论下, 检测不是必需的, 检测只是人为地设定阈值对目标是否存在进行判断。其过程可以描述为[5]:

在雷达等系统中, 检测常采用Neyman-Pearson准则, 即在系统满足一定虚警概率的情况下, 使得系统的检测概率最大。由此, 检测阈值γk要根据系统要求的虚警概率进行设定。

在检测到目标的情况下, 就需要对目标的状态进行估计, 在贝叶斯理论中, 常见的估计方法有2种: ①按照最大后验概率方法进行估计; ②按照贝叶斯最小均方误差方法进行估计。方法①对目标的状态进行估计, 即选择使得目标状态后验概率最大的目标状态作为估计值, 可以表示为:

按照方法②对目标状态进行估计, 即选择使得贝叶斯均方误差最小的目标状态作为目标状态的估计值, 可以表示为:

按照式 ( 9) 对目标状态进行估计需要进行积分或求和运算, 按照式 ( 8) 对目标状态进行估计则只需要确定目标状态后验概率密度函数最大值所对应的位置。

2 目标观测模型

从式 ( 4) 和式 ( 6) 中可以看出, 为了迭代计算目标状态的后验概率密度函数, 首先需要根据目标的动态方程对目标的状态进行预测, 接着需要根据目标的观测模型及实际观测值对预测进行修正。由此可见, 目标动态模型和观测模型直接影响着系统的性能。目标的动态方程常建立在直角坐标系下, 根据目标的机动能力按照匀速运动模型或匀加速运动模型建立动态方程[11,12]。为了使观测数据与目标的观测模型相吻合, 需要针对观测数据产生的机理建立目标的观测模型, 以避免由于目标观测模型失配所带来的损失。下面具体分析匹配滤波器的时频域特性对目标回波信号的影响。

在雷达等系统中, 为了充分利用时域信号的能量, 使输出某时刻的瞬时信噪比最大, 需要采用匹配滤波器[11]。信号通过匹配滤波器的输出不为冲激函数, 因此目标通过匹配滤波器的输出也会扩展到多个距离—多普勒单元, 利用目标在这些距离—多普勒单元上的观测, 以及它们位置的相关信息, 可以改善系统的检测和跟踪性能, 避免由于观测数据与观测模型不匹配所带到的性能损失。模糊函数表示匹配滤波器的输出, 对于带宽为5 MHz, 时宽为20μs的线性调频信号, 时延为0时的模糊函数如图1所示, 频延为0时的模糊函数如图2所示。

从图1和图2中可以看出, 目标经匹配滤波器后会分布在不同的距离—多普勒单元, 充分考虑它们之间的关系可以使系统的观测模型与实际观测数据更加吻合, 从而充分利用目标的散射能量, 提高系统的检测性能。

下面以目标在时域或距离维的扩展为例进行说明, 与目标在频域维的扩展情况类似。为了较充分利用目标在时域上分布的散射能量, 同时不显著增加系统的计算复杂度, 考虑到匹配滤波器的时域特性, 仅需考虑匹配滤波器输出主瓣的范围。考虑到信号的分辨能力, 在时域上常按照距离分辨力对匹配滤波后的信号进行采样。由此, 仅需考虑与目标位置最近的2个时域采样单元。在这2个时域采样单元上, 目标回波幅度受匹配滤波器的输出特性决定。设定匹配滤波器的时域输出归一化波形为so ( t) , 其对应的等效距离维输出波形为so ( R) , 目标位置为RT, 离目标位置最近的2个时域采样单元对应的距离分别为R1和R2, 有R1≤RT≤R2, 则目标在采样单元R1和R2对应的幅度调制因子分别为so ( RT- R1) 和so ( R2- RT) 。匹配滤波器输出为sinc函数时, 按照距离分辨力进行采样示意图如图3所示。

3 处理流程

匹配滤波器对波形相同而幅值不同的时延信号具有自适应性, 而对频移信号不具有自适应性[11]。结合匹配滤波器的上述性质及信号的模糊函数, 为了使得观测数据与目标的观测模型相吻合, 目标状态中的速度分量有2种处理方式。

第1种处理方式是直接结合匹配滤波器的时频域特性计算式 ( 6) 中的似然函数值p (zkx k) , 处理方式如图4所示。

图4中, 系统处理的对象为脉冲压缩后的信号, 针对匹配滤波器的时频域特性, 在不同的目标状态下计算观测似然函数。利用匹配滤波器的时频域特性和扫描方向图建立观测模型, 能够使观测数据与观测模型更吻合, 从而减小由于模型失配所带来的性能损失。

第2种处理方式利用目标状态中的速度分量调整匹配滤波器的冲激响应, 避免由于目标运动所带来的失配损失; 之后利用匹配滤波器的时域特性在不同的目标状态下计算似然函数, 如图5所示。

图5中, 系统处理的对象为脉冲压缩前的信号, 在不同的目标状态下, 利用目标状态的速度分量调制匹配滤波器的冲激响应, 使匹配滤波器与目标回波信号匹配, 从而避免由于匹配滤波器失配所带来的损失。脉冲压缩之后, 针对匹配滤波器的时域特性计算观测似然函数。

第1种处理方式直接处理脉冲压缩后的数据, 根据匹配滤波器的时频域特性建立目标的观测模型, 使观测数据能够与观测模型很好地吻合, 减小系统性能的损失。

第2种处理方式需要处理脉冲压缩前的数据, 根据目标的状态设置匹配滤波器的参数, 使匹配滤波器能够与回波信号匹配, 减小由于回波信号频移所带来的失配损失, 提高回波信号的信噪比; 另外, 对于时域频域耦合的信号 ( 如线性调频信号) , 能够在脉冲压缩阶段除去时域和频域之间的耦合性, 使距离测量更加准确。

两种处理方式都针对回波信号的模型建立目标的观测模型, 使观测数据能够与目标的观测模型较好地吻合, 减小由于模型失配所带来的损失。第2种处理方式相对于第1种处理方式, 能够减小由于回波信号频移所带来的匹配损失, 进一步提高脉冲压缩后的信噪比。另外, 第2种处理方式在脉压阶段也解决了回波信号时域和频域耦合的问题; 而对于第1种处理方式, 则需要在计算观测似然函数阶段考虑回波信号时域频域的耦合。然而, 第2种处理方式需要在目标所有的状态情况下对原始观测数据进行脉冲压缩处理, 计算量较大。

4 数值仿真及分析

下面对上述模型进行仿真分析, 以验证系统的跟踪性能。在贝叶斯检测前跟踪算法中, 从式 ( 2) 中可以看出, 观测与目标状态之间呈现高度非线性关系, 由此很难用一般的跟踪滤波方法进行处理。对于贝叶斯检测前跟踪算法, 为了求解式 ( 4) 和式 ( 6) , 通常有两种近似方法: 第1种方法是采用网格近似方法, 即将目标的状态离散化, 由此式 ( 4) 所示的积分运算就可以转变为求和运算; 第2种方法是采用粒子滤波方法, 用大量粒子的随机运动来近似目标状态的后验概率密度[4]。

下面采用粒子滤波的方法对目标状态的后验概率密度进行近似, 在SIR滤波器的基础上, 采用MCMC方法解决SIR粒子滤波器中粒子的多样性, 采用Auxiliary SIR滤波器的方法使得粒子更集中在目标可能出现的位置附近[4]。

考虑距离维和多普勒维, 目标在雷达的射线范围做直线运动, 速度分量受到一个随机扰动, 目标在连续时间状态下的状态方程可以表示为:

式中, 表示距离和多普勒;w ( t) 为零均值的高斯系统噪声。设定目标的初始状态为[100 km, 250 m/s]', w ( t) 的功率谱密度设为5。

发射带宽为5 MHz, 时宽为20μs的线性调频信号, 脉冲压缩前按照2倍信号带宽进行采样, 脉冲压缩后按照1倍带宽进行重采样, 由此可知, 距离分辨率为30 m。

脉冲压缩后的信噪比设定为SNR =9 dB, 考虑匹配滤波器时频域特性前后, 贝叶斯检测前跟踪算法目标位置的均方根误差如图6所示, 其中实线表示考虑匹配滤波器时频域特性时目标位置的均方根误差, 虚线表示不考虑匹配滤波器时频域特性时目标位置的均方根误差。

从图6中可以看出, 考虑匹配滤波器时频域特性时目标位置的均方根误差要小于不考虑匹配滤波器时频域特性时目标位置的均方根误差, 这主要是因为考虑匹配滤波器时频域特性时, 目标的观测模型与处理的数据更加吻合, 减小了由于目标观测模型失配引起的性能损失, 改善了贝叶斯检测前跟踪算法的性能。

5 结束语

目标通过匹配滤波器的输出会扩展到多个距离—多普勒单元, 利用目标在这些距离—多普勒单元上的观测以及它们位置的相关信息, 可以改善贝叶斯检测前跟踪算法的性能。本文详细分析了观测信号的产生机理, 根据匹配滤波器时频域特性建立合适的观测模型, 给出了利用匹配滤波器时频域特性的贝叶斯检测前跟踪算法的处理流程, 从而避免由于观测数据与观测模型不匹配所带来的性能损失, 改善了贝叶斯检测前跟踪算法的性能。本文的结论对贝叶斯检测前跟踪算法在实际工程中的应用有一定的指导意义。

参考文献

[1]CARLSONB D, WILSON S L.Search Radar Detection and Track with the Hough Transform Part III:Detection Performance with Binary Integration[J].IEEE Trans on AES, 1994, 30 (1) :116-125.

[2]李涛, 吴嗣亮, 曾海彬, 等.基于动态规划的雷达检测前跟踪新算法[J].电子学报, 2008, 36 (9) :1 824-1 828.

[3]TONISSEN S, SHALOM Y B.Maximum Likelihood TrackBefore-Detect with Fluctuating Target Amplitude[J].IEEE Trans on AES, 1998, 34 (3) :796-809.

[4]RISTIC B, GORDON N.Beyond the Kalman Filter:Particle Filters for Tracking Applications[M].London:Artech House, 2004.

[5]吴孙勇, 廖桂生, 杨志伟.改进粒子滤波的弱目标跟踪与检测[J].宇航学报, 2010, 31 (10) :2 395-2 401.

[6]MORELANDE M R, KASTELLA K.A Bayesian Approach to Multiple Target Detection and Tracking[J].IEEE Trans on SP, 2007, 55 (5) :1 589-1 604.

[7]ORLANDO D, LOPS M, RICCI G.Track-Before-Detect Strategies for STAP Radars[J].IEEE Trans on SP, 2010, 58 (2) :933-938.

[8]BUZZIS, VENTURINO L.Track-Before-Detect Procedures for Early Detection of Moving Target from Airborne Radars[J].IEEE Trans on AES, 2005, 41 (3) :937-954.

[9]焦子朋, 耿哲, 周思玥.MTI在连续波雷达中的应用研究[J].无线电工程, 2012, 42 (11) :62-64.

[10]丁静.稳健数字波束形成算法研究[J].无线电通信技术, 2012, 38 (6) :30-32.

[11]吴顺君, 梅晓春.雷达信号处理与数据处理技术[M].北京:电子工业出版社, 2008.

强跟踪滤波滤波器 篇6

传统的火灾检测技术一般都是只基于概率模型, 即直接对火灾疑似区域进行各个特征的检测来综合判断火灾的发生。这种检测技术虽然实现简单,但由于环境多变性的影响,误检率较高。对于火灾监控系统而言,由于需 24 小时不间断运行,环境特别是光线的干扰对算法检测效果会产生较严重影响。此外,实时检测系统对每帧图像的处理时间有一定限制,因此,算法复杂度不能太高,但算法简单易导致较高的误报率。针对准确性与实时性之间的矛盾,本文先对疑似火灾区域进行Kalman跟踪,由于光线的变化较为缓慢,必然导致Kalman跟踪在时间上的不连续,而火焰一旦产生就会持续燃烧,Kalman跟踪在时间上是连续的。依据火焰燃烧的连续性, 并对疑似火灾区域进行火焰特征检测, 进而经过连续多帧的综合判断来识别火灾。实验结果表明,本文提出的算法响应时间短,抗干扰能力强,能够满足实际火灾监测的要求。

1 Kalman 滤波器

Kalman滤波器最初于1960年提出,其基本思想是根据系统的历史测量值,建立最大化这些历史测量值的后验概率的系统状态模型。Kalman滤波以目标的位置、速度和加速度来描述目标状态矢量,根据递推的方法,不断更新目标状态来预测目标。

Kalman滤波器假设被建模的系统是线性的[1], “线性”的意思是指Kalman滤波器定义中的不同阶段可以用矩阵乘法来表示,这决定了Kalman滤波器跟踪快速移动位移量较大的物体的效果不是很好,但火焰的移动位移较小,适合用Kalman 滤波器进行跟踪。

Kalman滤波步骤[2,3]如下:

Step1 预测。设输入的n维目标状态矢量为x,则计算先验估计xk的公式为:

$x{}_{\overline{k}}=Fx{}_{k-1}+\mathit{B}u{}_k+w{}_k(1)$

从式(1)中可以看出,将时间 k 描述为时间 k-1 的函数。其中F是一个n×n矩阵,也被称为传递矩阵。uk的作用是允许外部控制此系统,由表示输入控制的c维向量组成。B是一个联系输入控制与状态改变的n×c矩阵。wk被称为过程噪声,是一个影响状态系统的随机外部变量。

Step2 计算测量值zk。公式如下:

$z{}_k=\mathit{{\rm H}}{}_{k}x{}_{\overline{k}}+v{}_k(2)$

这里,Hk为 m×n 维的测量矩阵,vk为测量误差。从式(2)中可以看到,在时刻 k获得的新的测量值是与旧时刻 k-1 向前映射到时刻 k 的模型融合得到的, 而不是直接与旧时刻 k-1 的模型融合。

Step3 更新误差协方差矩阵以及计算Kalman更新率误差协方差矩阵在时间k的先验估计由其在时刻k-1的值来更新。

$\mathit{{\rm P}}{}_{\overline{k}}=\mathit{F}\mathit{{\rm P}}{}_{k-1}\mathit{F}{}^{{\rm T}}+\mathit{Q}{}_{k-1}(3)$

其中,F仍为传递矩阵,FT为其转置矩阵,Qk-1为过程噪声wk的协方差矩阵。

Kalman 更新率计算公式为:

$\mathit{{\rm K}}{}_k=\mathit{{\rm P}}{}_{\overline{k}}\mathit{{\rm H}}{}_k^t(\mathit{{\rm H}}{}_k\mathit{{\rm P}}{}_{\overline{k}}\mathit{{\rm H}}{}_k^t+\mathit{R}{}_k){}^{-1}(4)$

其中,${\rm P}{}_{\overline{k}}$为更新过的误差协方差,Hk仍为测量矩阵,Rk为测量误差 vk的协方差矩阵。

Step4 计算修正后的目标状态矢量和Pk最优值 修正后的目标状态矢量由下式得出:

$x{}_k=x{}_{\overline{k}}+\mathit{{\rm K}}{}_k(z{}_k-\mathit{{\rm H}}{}_{k}x{}_{\overline{k}})(5)$

从式(5)中可以看出,最终的估计值是采用上一帧的估计值和当前帧的观测值综合计算得到的,这种实时递推过程充分利用了以前的状态信息,故估计值准确性较高。

求得更新率Pk最优值,以便进行下一次迭代:

$\mathit{{\rm P}}{}_k=({\rm I}-\mathit{{\rm K}}{}_k\mathit{{\rm H}}{}_k)\mathit{{\rm P}}{}_{\overline{k}}(6)$

Step5 以xk为起点,进行下一次预测。

2 基于Kalman 滤波器运动跟踪算法的火灾识别

基于Kalman滤波器运动跟踪算法的火灾识别方法利用 Kalman 滤波器跟踪运动区域,从而可以知道被跟踪区域是否连续存在, 并对连续多帧存在的运动区域进行特征检测来判别火焰,从而有效地排除干扰。具体算法流程如图1所示。

2.1 高斯前景检测

运动目标检测就是将运动目标从视频图像中提取出来, 本文通过运动目标检测方法把火焰区域从视频图像中检测出来。 目前常用的运动目标检测算法有帧间差分法、 背景差分法、光流法和高斯背景建模法。差分法、背景差分法虽然相对简单,但存在着目标轮廓检测不完整和目标相关点保留较少的问题; 而光流法更适合摄像机运动的情况, 通过求解偏微分方程求得图像序列的光流场,从而预测摄像机的运动状态。但是对于摄像机固定的情形,由于光流法的复杂性, 往往难以实时的计算, 所以本文采用高斯背景建模的方法。 通过对背景建模,对一幅给定图像分离前景和背景,一般来说,前景就是运动物体,从而达到对运动物体检测的目的。

高斯建模法首先选取N帧原始图像作为获取背景帧图像的来源的训练图像。高斯模型认为, 对于一个背景图像, 其特定像素亮度的分布满足高斯分布, 即背景图像B(x,y)点的亮度满足:

I(x,y)～N(u,σ) (7)

其中,u和σ表示每个背景像素的平均值和方差,且每一点的高斯分布是独立的。将背景中的每个像素按照高斯分布模型建模, 通过N帧原始图像的训练获得其参数并不断更新其分布参数, 据此更新背景图像来获取初始背景帧[4]。

高斯背景建模步骤[5]如下:

Step1 初始化背景帧图像。选取第一帧图像每个像素点的灰度值作为均值u,方差σ设为0。

Step2 更新背景模型。对于一幅输入图像G,若某点的灰度值G(x,y)满足式(8),则认为(x,y)是背景点,反之是前景点。

$Exp{}^{\frac{-(G(x,y)-u(x,y)){}^2}{2x\sigma {}^2}}＞{\rm T}(8)$

Step3 更新每个像素点的参数。随着时间的变化,背景图像也会发生缓慢的变化,这时要不断更新每个像素点的参数。

u(t+1,x,y)=a×u(t,x,y)+(1-a)×I(x,y) (9)

式中,a为更新参数,表示背景变化的速度,一般情况下,不更新σ(实验中发现更不更新σ,效果变化不大) 。

2.2 轮廓标记

经过高斯背景建模后,得到了二值化的前景图像, 并在此前景图像上提取运动区域的轮廓,然后给各个轮廓添加标记。得到标记的轮廓后,计算每个轮廓的中心点坐标以及轮廓的宽和高等信息,以便作为跟踪过程的参考状态。

2.3 Kalman 跟踪运动区域

经过标记的轮廓便于利用 Kalman滤波器进行跟踪。 设初始输入状态向量为(x,y,w,h,0,0),其中x和y表示运动区域轮廓的中心点坐标,w 和h表示轮廓的宽和高。

传递矩阵为:

$\left[\begin{array}{cccccc}1& 0& 0& 0& 1& 0\\ 0& 1& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right](10)$

测量矩阵为:

$\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right](11)$

协方差阵对滤波的影响不大,可以用对角阵作为它的初始值[6]。

过程噪声的协方差矩阵为:

$\left[\begin{array}{cccccc}10{}^{-5}& 0& 0& 0& 0& 0\\ 0& 10{}^{-5}& 0& 0& 0& 0\\ 0& 0& 10{}^{-5}& 0& 0& 0\\ 0& 0& 0& 10{}^{-5}& 0& 0\\ 0& 0& 0& 0& 10{}^{-5}& 0\\ 0& 0& 0& 0& 0& 10{}^{-5}\end{array}\right](12)$

测量误差的协方差矩阵为:

$\left[\begin{array}{cccc}0.1& 0& 0& 0\\ 0& 0.1& 0& 0\\ 0& 0& 0.1& 0\\ 0& 0& 0& 0.1\end{array}\right](13)$

初始误差协方差矩阵为:

$\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right](14)$

设外部输入控制uk=0,即此系统不受外部控制,自动跟踪轮廓目标区域,把初始矩阵代入式(1)-(6)中进行跟踪。

2.4 火焰特征检测

根据火灾初期火焰表现出来的特征信息、 特征提取的可行性与有效性, 本文选取了火焰颜色特征、圆形度、火焰运动特征来逐级识别火焰。

①火焰颜色判断

在彩色图像中, 火焰的内核呈现亮白色, 向外随着温度的降低颜色逐渐由黄转红。 因此,火焰的红色分量和亮度是区别火焰的一个显著特征。 同时, 火焰颜色分量之间存在内在关系,红色分量大于等于绿色分量,绿色分量又大于等于蓝色分量[7];另外,火焰区域图像的饱和度较高,亮度值大。经过多次测试比较:本文选取下列颜色判定公式来判定符合火焰颜色的特征点:

R(x,y)≥105 (15)

R(x,y)≥G(x,y)≥B(x,y) (16)

S(x,y)≥0.2 (17)

$\begin{array}{l}S(x,y)=\frac{\max -\min }{\max }(18)\\ {\rm I}(x,y)=\frac{R(x,y)+G(x,y)+B(x,y)}{3}(19)\\ \max =\max (R(x,y),G(x,y),B(x,y))(20)\\ \min =\min (R(x,y),G(x,y),B(x,y))(21)\end{array}$

上述公式中,R(x,y),G(x,y),B(x,y),S(x,y),I(x,y)分别代表图像坐标为(x,y)处像素点的红绿蓝颜色分量值以及饱和度与亮度值。扫描原始图像运动区域内的每一个像素点, 若该像素点同时满足上述公式, 则判定该像素点为火焰疑似点。 统计该运动区域内火焰疑似点的总数,计算火焰疑似点占该运动区域总像素点总数的百分比。经过大量统计,若此百分比大于0.1,可将该区域记为疑似火焰区域,进行下一步判断。

②圆形度

对经过颜色判断的疑似火焰区域进行圆形度的检测。

圆形度是反映物体或区域的形状复杂度的特征量。圆形度的定义为:

$e=\frac{4\times {\rm P}{\rm I}\times S}{L{}^2}(22)$

式中,S为物体或区域的面积,L指物体或区域的周长。

周长是指区域边界的长度。 本文中所求的周长是指运动轮廓边界白点的个数, 通过统计区域的白点总数可求得面积。不规则物体的圆形度一定小于1。在经过大量实验后, 假定火焰圆形度在0.2～0.6的区间内。通过遍历整幅图像,计算出每个轮廓的圆形度,如果圆形度值在范围内,则判定这个区域为疑似火焰的区域。

圆形度作为早期表征火焰特性的判据是十分有效的, 可以排除一些规则发光体 (如白炽灯、手电筒等)的干扰,从而为下面的判据减小计算量。

③火焰运动特征

由于火焰燃烧过程中形状的不规则性, 使得火焰运动区域的中心在水平和垂直方向上来回抖动。而人的位移则在水平和垂直方向上能够保持很好的连续性。因此。通过记录连续N帧图像每个运动区域的中心点坐标位移在水平和垂直方向上的偏移量, 可识别出火焰运动区域并排除掉呈现规则运动方向的人的干扰。轮廓中心(x,y)计算公式为:

$X=\frac{{\displaystyle \sum _{i=1}^{{\rm M}\times {\rm N}}x}{}_i}{{\rm M}\times {\rm N}},Y=\frac{{\displaystyle \sum _{i=1}^{{\rm M}\times {\rm N}}y}{}_i}{{\rm M}\times {\rm N}}(23)$

其中,M和N分别表示轮廓的宽和高。

用(Xn,Yn)、( Xn-1,Yn-1)分别表示相邻两帧图像同一轮廓的中心点坐标, (DXn,DYn)表示相邻两帧图像同一轮廓中心点坐标的差值, 设置计数器CT1和计数器CT2, 若 DXn则 CT1加1, 若DYn则 CT2加 1,取N=10,即 :

$\{\begin{array}{l}DX{}_n=X{}_n-X{}_{n-1}\\ DY{}_n=Y{}_n-Y{}_{n-1}\\ \text{i}\text{f}(DX{}_n＞0)C{\rm T}{}_1++\\ \text{i}\text{f}(DY{}_n＞0)C{\rm T}{}_2++\end{array}n=2,3,\cdots ,10(24)$

由于火焰抖动导致其形状的不断变化,故DXn和DYn总在大于零和小于零之间来回摆动。而人的位移总是朝着一个方向运动,故在连续10帧图像中,其DXn和DYn总有一个橫大于零或小于零。经过大量实验,若CT1和CT2同时满足式:

(2<CT1<8)&&(2<CT2<8) (25)

可认为该区域满足火焰特征。

3 实验分析

为了验证本文算法的有效性,本文参考了 ISO/WD 7240-29《火灾探测报警系统 第29部分:图像型火灾探测器》进行火灾检测试验。测试环境如下:

处理器:Intel(R) Core(TM) i5-2320 CPU @3.00GHz。

内存:8G。

操作系统:Windows 7 旗舰版,64 位操作系统。

摄像头:540 线,距离火源 10m。

录像机:4 路 D1 格式,H.264 硬盘录像机。

实验场所:室内。

温度:27℃。

为了验证本文算法能有效排除不同条件的干扰, 实验环境中加入了光线的明暗变化以及人的干扰。

由于录像机传输过来图像为 YUV 格式,故需经过颜色空间转换为 RGB格式,以便进行后续的火焰颜色判断。转换公式如下所示:

R=1.164×(Y-16)+1.596×(V-128)

G=1.164×(Y-16)-0.392×(U-128)-0.813×(V-128)

B=1.164×(Y-16)+2.017×(U-128) (26)

对转换之后的图像进行后续处理,所得实验结果如表1所示。

从图 2(a)图中可以看出,实验环境中存在树影的晃动和光线的明暗变化,图中的椭圆区域为检测出的火焰区域。图2(b)中,在夜晚光线较暗的环境下检测火焰,且环境中存在人的干扰。

大量的实验结果表明,本文提出的算法能正确地识别火焰,且响应时间短,能基本排除环境干扰,实现准确性与实时性的统一。

4 结束语

本文通过把 Kalman 滤波器的运动跟踪算法应用于火焰检测,实现了准确性与实时性相统一的火焰识别系统。实验结果表明,通过跟踪运动区域,并结合颜色检测、圆形度检测、火焰运动特征分析等步骤, 本文提出的算法可以较好地抑制环境中光线变化的干扰, 降低实时监控系统的误报率,且能够准确检测出火焰,具有很强的实用性。

参考文献

[1]Khemiri k.Robust fault and state estimation for linear discrete-timesystems with unknown disturbances using PI Three-Stage KalmanFilter[J].Communications,Computing and Control Applications,2011,3:1-7.

[2]Chen S Y.Kalman Filter for Robot Vision:A Survey[J].Industrialelectronics,2012,11(11):4409-4420.

[3]虞旦,韦巍,张远辉.一种基于卡尔曼预测的动态目标跟踪算法研究[J].光电工程,2009,36(1):52-57.

[4]刘亚,艾海舟,徐光佑.一种基于背景模型的运动目标检测与跟踪算法[J].信息与控制,2002,31(4):315-319.

[5]康晓晶,吴谨.基于高斯背景建模的目标检测技术[J].液晶与显示,2010,25(3):454-459.

[6]Facchini M.Distributed optical fiber sensors based Brillouin scatter-ing[D].Lausanne,Switzerland:Swiss Federal Institute of Technolo-gy,2001.

强跟踪滤波滤波器 篇7

1 视频跟踪问题描述

最优贝叶斯估计的基础由式 (1) 和式 (2) 组成, 式 (1) 为预测方程、式 (2) 为更新方程。因为有高维积分运算的存在, 因而式 (1) 和式 (2) 很难求出解析解, 通常情况下都是通过Kalman滤波方法、基于网格的滤波方法、Unscented Kalman滤波方法及Monte Carlo方法等近似求解。

2 Kalman滤波方法

Kalman滤波方法主要是将一些特征量, 如目标的形状、速度、颜色、纹理等当作状态量, 且对这些状态量进行分析, 其一般过程是首先估算某一时刻的状态, 然后将估算的状态进行反馈, 具体反馈是通过测量变量实现的。

目标的状态方程为式 (4) 。

目标的观测方程为式 (5) 。

式 (4) 与式 (5) 中的Xk、Zk分别为时刻的状态向量和观测向量, A与H分别为状态转移矩阵和观测矩阵, Wk-1、Vk分别代表k-1和k时刻独立的且没有相关性的状态噪声和观测噪声, 且都满足均值为零的高斯分布, 其协方差矩阵为式 (6) 和式 (7) 。

滤波使各个时刻目标状态向量的后验估计值误差协方差最小, 主要包括预测和修正两个部分, 其中, 预测部分是对先验估计误差协方差和状态进行预测;修正部分是对滤波增益的计算及使用滤波增益修正后验估计误差协方差与状态预测值, 通过运算得到最小化的后验估计误差协方差与状态向量的后验估计, 预测部分与修正部分的表达式如下所示。

状态预测方程 (先验估计) 为式 (8) 。

先验估计误差协方差为式 (9) 。

其中, 式 (9) 中Q为状态噪声协方差矩阵。

经观测值修正过的后验估计为式 (10) 。

后验估计误差的协方差为式 (11) 。

3 计算结果与分析

图1为Kalman视频跟踪流程图。图2为视频中的五帧图像和相应的跟踪结果。由图2可看出, Kalman滤波算法可预测运动目标下一帧的位置, 实现对运动目标的跟踪。

4 结语

针对视频目标跟踪, 阐述基于Kalman滤波的跟踪算法, 并进行仿真分析, 结果实现对运动目标的跟踪。

参考文献

[1]Peterfreund N.Robust Tracking of Position and Velocity with Kalman Snakes[J].IEEE Computer Society, 2002 (6) :564-569.

[2]Nikos Paragios, Rachid Deriche.Geodesic Active Contours and Level Sets for the Detection and Tracking of Moving Objects[J].IEEE Trans Pattern Analysis and Machine Intelligence, 2000 (3) :266-280.

[3]Michael Isard, Andrew Blake.Contour Tracking by Stochastic Propagation of Conditional Density[J].European Conference on Computer Vision Cambridge, 1996 (4) :343-356.

[4]Nummiaro K, Koller-Meier E, Van Gool L.An Aadaptive Color-based Particle Filter[J].Image and Vision Computing, 2003 (1) :99-110.

[5]Vignon D, Andrews L, Robert J.General Purpose Real-time Object Tracking Using Hausdorff Transforms[M].France:ESIA, 2002:1-6.

[6]Polana R, Nelson R.Low Level Recognition of Human Motion (or How to Get Your Man without Finding His Body Parts) [J].IEEE Xplore, 1999 (10) :77-82.

[7]Jang D-S, Choi H-I.Active Models for Tracking Moving Objects[J].Pattern Recognition, 2000 (7) :1135-1146.

基于粒子滤波的声源方位跟踪算法 篇8

基于麦克风阵列的声源定位与跟踪技术可广泛用于视频会议[1]、语音增强[2]、智能机器人[3] 等领域。例如在视频会议系统中,基于麦克风阵列的声源定位技术可为摄像头提供说话人的方位,使摄像头自动对准说话人。当说话人走动时,还可使用跟踪算法对其进行跟踪。粒子滤波PF(Particle Filtering)是在贝叶斯滤波框架下基于Monte Carlo 采样的统计滤波方法,可解决非线性、非高斯的滤波问题,声源跟踪就属于这类问题。近年来,基于粒子滤波的声源跟踪算法已成为研究的热点。

文献[4]采用Langevin方程来模拟声源的运动,首次将粒子滤波用于解决噪声和混响环境中的声源跟踪问题。文献[5]在此基础上做了改进,并总结了几种声源跟踪的方法。此后,为了提高跟踪精度,国内外的研究人员又提出了多种改进算法[6,7,8]。跟踪算法中的似然函数对跟踪性能起着决定性的作用,声源定位算法是构造似然函数的基础,对跟踪性能起着关键性的作用。相位变换加权的可控响应功率SRP-PHAT(Steered Response Power-Phase Transform)声源定位算法[9]在混响环境中有较强的鲁棒性,定位精度较高,在跟踪算法中常用作定位函数。文献[4]提出的声源运动模型理论上可用于三维空间中的声源跟踪,但为了达到较高的跟踪精度,需要采用分布式的大型麦克风阵列。此外三维跟踪的计算量很大,实现实时跟踪比较困难。

在麦克风阵列的很多实际应用中,只需要估计声源的方位,当声源移动时,也只需要实现声源的方位跟踪。如用于视频会议和语音增强的麦克风阵列,通常只要给出说话人的方位。本文提出一种基于粒子滤波的声源方位跟踪算法。该算法在Langevin方程的基础上构建声源的方位动态模型,采用远场条件下的SRP-PHAT函数构造似然函数,有效地实现了真实环境中说话人方位的跟踪。

1 贝叶斯滤波

声源跟踪问题通常指在混响环境中根据M个麦克风的接收信号,实时估计出当前声源的位置或方位。接收信号按时间顺序分成连续的若干帧,帧长为L,则第m个麦克风(m=1,2,…,M)的第t帧信号为:

其中“T”表示转置。 M个麦克风的第t帧信号组成一个M×L的矩阵

在球坐标系统下,声源位置可表示为q=[r,θ,ϕ]T,其中r,θ,ϕ分别表示距离,水平角和仰角。假设目标处在阵列远场,则定位算法只能估计声源的方位,也就是声源的到达方向角DOA(Direction of Arrival),用l=[θ,ϕ]T来表示。声源方位跟踪就是估计声源在当前时刻的l。在远场假设下,目标在第t帧时的状态矢量可表示为:

$\alpha {}_t=[\theta {}_t,\text{ϕ}{}_t,\overline{\theta }{}_t,\overline{\varphi }{}_t]{}^{\text{Τ}}$ (3)

其中[θt,ϕt]T表示声源方位,[$\overline{\theta }$t,$\overline{\varphi }$t]T表示声源的角速度。给定状态矢量α,我们用lα表示相应的声源方位。Yt表示目标状态的观测值,该观测值由接收信号Xt,通过定位函数转换得到,即:

其中f(·)为定位函数,与所采用的定位方法有关。假设状态转移是一个Markov过程,则声源跟踪问题可以用如下两个方程来描述:

Yt=h(αt,vt) (6)

其中g(·)和 h(·)可能是非线性函数,ut和vt可能为非高斯的噪声矢量。用Y1:t={Y1,…,Yt}表示从初始时刻直到t时刻的观测值,贝叶斯滤波问题就是在每个时刻t,利用所获得的观测值Y1:t估计状态αt的后验概率密度函数p(αt|Y1:t),通过求数学期望,得到t时刻的声源状态估计$\widehat{\alpha }$t。贝叶斯滤波问题的解包括预测和更新两个步骤。假设t-1时刻的后验概率密度函数p(αt-1|Y1:t-1)是已知的,则t时刻的后验概率密度函数p(αt|Y1:t)可通过以下两个方程得到:

$p(\alpha {}_t|Y{}_{1:t-1})={\displaystyle \int \begin{array}{c}p(\alpha {}_t|\alpha {}_{t-1})\end{array}}p(\alpha {}_{t-1}|Y{}_{1:t-1})\text{d}\alpha {}_{t-1}$ (7)

其中p(αt|Y1:t-1)是状态一步预测概率密度函数,p(αt|αt-1)为状态转移概率密度,p(Yt|αt)为似然函数。

2 基于粒子滤波的声源方位跟踪

通常,方程式(7)和式(8)没有闭式解存在,除非方程式(5)和式(6)中的g(·)和 h(·)是线性函数, ut和vt是高斯噪声矢量,此时卡尔曼滤波是最优解。声源跟踪问题中,这些函数是非线性的,因此不能使用卡尔曼滤波求解,而粒子滤波则是解决这类非线性、非高斯问题的强有力的工具。粒子滤波是一种按Monte Carlo仿真实现递推贝叶斯滤波的技术,其关键思想是根据一组带有相应权重的随机样本来表示后验概率密度函数,而且基于这些样本和权重来计算估计值。状态空间的样本α${}_t^{(i)}$称为粒子,其相应权重为w${}_t^{(i)}$,i∈{1,…,N}。当样本数N非常大时,这些样本和权重就趋于真实地描述了后验概率密度函数,因而粒子滤波就成为了跟踪问题的最优贝叶斯解[10]。

方程式(5)的具体实现称为声源动态模型,文献[4]提出用Langevin 方程来模拟声源的运动,该模型假设,在直角坐标系中,声源在每个坐标轴方向的随机运动是独立同分布的。在二维平面运动声源跟踪实验中,使用Langevin模型得到了较高的跟踪精度。远场条件下,声源运动过程中,水平角和仰角的变化也可以认为是独立同分布的随机过程。譬如声源水平角的变化可描述为:

$\dot{\theta }{}_t=a{}_{\theta }\dot{\theta }{}_{t-1}+b{}_{\theta }F{}_{\theta }$ (9)

$\theta {}_t=\theta {}_{t-1}+\Delta {\rm T}\dot{\theta }{}_t$ (10)

aθ=exp(-βθΔT) (11)

$b{}_{\theta }=\overline{v}{}_{\theta }\sqrt{1-a{}_{\theta }^2}$ (12)

其中,Fθ是均值为0,方差为1的正态分布随机变量,$\Delta {\rm T}=\frac{L}{f{}_s}$,是两次声源方位估计的时间间隔(L是信号帧长,fs是信号采样频率),$\overline{v}$θ是稳态角速度,βθ为常数。参照文献[4]的建议,设置βθ=10Hz,仰角方向的对应参数βφ=10Hz。说话人在运动过程中,水平角变化较快,而仰角仅当声源高度变化(如说话人弯腰)时才有明显变化。根据此特点,设置$\overline{v}{}_{\theta }=1\text{r}\text{a}\text{d}/\text{s}$,$\overline{v}{}_{\text{ϕ}}=0.4\text{r}\text{a}\text{d}/\text{s}$。

定位函数将接收到的原始信号转变成目标状态观测值,对跟踪性能起着非常重要的作用。如前文所述,SRP-PHAT定位精度较高,跟踪算法中常用作定位函数。给定声源方位l,SRP-PHAT的函数表达式为:

$f(l‚X{}_t)={\displaystyle \sum _{l=1}^{{\rm M}}}{\displaystyle \sum _{m=l+1}^{{\rm M}}\widehat{R}}{}_{l\text{m}}[\text{τ}{}_{lm}(\text{l})]$ (13)

其中$\widehat{R}{}_{lm}[\tau {}_{lm}(l)]$为第l个和第m个麦克风的接收信号的广义互相关函数(GCC: Generalized Cross-Correlation),其表达式为:

$\widehat{R}{}_{lm}(\tau )=\frac{1}{{\rm K}}{\displaystyle \sum _{k=0}^{{\rm K}-1}\frac{X{}_l(k)X{}_m^{*}(k)}{\left|X{}_l(k)X{}_m^{*}(k)\right|}}e{}^{j\omega \tau }$ (14)

其中Xm(k)是xm(t)的FFT,K为FFT点数,ω是模拟角频率,“*”表示取共轭,τlm(l)为假想声源到第l个和第m个麦克风的到达时间差TDOA(Time Difference of Arrival)。在远场假设下,麦克风阵列接收信号为平面波,τlm(l)可用下式来计算:

$\tau {}_{lm}(l)=\frac{\zeta {}^{{\rm T}}(r{}_l-r{}_m)}{c}$ (15)

其中ζ为声源的单位方向矢量,其表达式为:

ζ=[cosϕcosθ, cosϕsinθ, sinϕ]T (16)

rm=[x y z]T为第m个麦克风在直角坐标系中的坐标矢量,c为空气中的声速(约为342m/s)。

似然函数的作用是评价粒子权重,定位函数是声源方位的连续函数,因而可用作似然函数,我们把它称为伪似然PL(Pseudo Likelihood)函数。我们采用的伪似然函数为:

这里取Yt(lα)和0的最大值的作用是为了保证似然函数非负,η∈R+,其作用是使似然函数更尖锐,从而更适合于声源跟踪[5]。将式(13)代入式(4),再将该式代入式(17)即可构造出对应SRP-PHAT的似然函数。基于粒子滤波的声源方位跟踪算法流程:

产生一组初始粒子{α${}_0^{(i)}$,i=1:N},赋予均匀权重{w${}_0^{(i)}$=1/N,i=1:N},对接收到的每一帧新的数据执行以下步骤:

① 根据状态转移方程预测新的粒子组,$\tilde{\alpha }{}_t^{(i)}=g(\alpha {}_{t-1}^{(i)},u{}_t)$;② 通过定位函数将原始数据转变成目标状态观测值,Yt=f(l,Xt);

③ 形成似然函数,$p(Y{}_t|\tilde{\alpha }{}_t)=F(Y{}_t,\tilde{\alpha }{}_t)$;

④ 根据似然函数计算新的粒子权重,$\tilde{w}{}_t^{(i)}=p(Y{}_t|\tilde{\alpha }{}_t^{(i)})$,然后将权重归一化,$w{}_t^{(i)}=\tilde{w}{}_t^{(i)}\cdot ({\displaystyle \sum _{i=1}^{{\rm N}}\tilde{w}}{}_t^{(i)}){}^{-1}$;

⑤ 计算当前时刻声源方位估计,$\widehat{l}{}_s(t)={\displaystyle \sum _{i=1}^{{\rm N}}w}{}_t^{(i)}l{}_{\tilde{\alpha }}^{(i)}$;

⑥ 按照权重,从现有粒子组$\{\tilde{\alpha }{}_t^{(i)},i=1$:N}中重采样,得到重采样粒子组{α${}_t^{(i)}$,i=1:N},然后将所有粒子置为均匀权重,w${}_t^{(i)}$=1/N;

⑦ 存储重采样后的粒子和它们的权重{α${}_t^{(i)}$,w${}_t^{(i)}$,i=1:N},回到①。

3 实验结果与分析

为了验证本文提出算法的有效性,我们使用取自瑞士IDIAP研究所的真实数据[11]来做声源跟踪实验。IDIAP提供了单声源、多声源、静态和动态等多种情景的录音。我们取其中编号为“seq11-1p-0100”的一组数据。该组数据是在一间普通小型会议室中录制的,麦克风阵列如图1所示。

该阵列为均匀圆阵,半径为0.1m,图中黑点表示麦克风,编号1-8。麦克风阵列放置于会议桌上,录制过程中,一个说话人面对阵列在会议桌旁边走边说,期间有弯腰动作,因此声源高度有变化。录制场景如图2所示,图中虚线表示说话人在地面走动的轨迹。

声源信号是一段男声英语语音,时长约6.6s,采样频率16kHz。IDIAP提供了说话人走动过程中嘴部的三维坐标,每40ms给出一组数据,根据这些数据可以得到声源的真实方位。为了便于与真实方位比较,跟踪算法中取信号帧长L=640点(40ms),帧之间不重叠,加Hanning窗,接收信号总共分为165帧。粒子数N=50,式(17)中的η=3。跟踪精度用均方根误差RMSE(Root Mean Square Error)来评价,定义为:

${\rm E}{}_{R{\rm M}SE}=\sqrt{\frac{1}{{\rm K}{}_s}{\displaystyle \sum _{t=1}^{{\rm K}{}_s}l{}_s(t)-\widehat{l}{}_{\text{s}}(t)}{}^2}$ (18)

其中,ls(t)表示第t帧时刻声源的真实方位,$\widehat{{\displaystyle l}}{}_s(t)$为其估计值,‖·‖表示求该矢量的2-范数, Κs为信号帧数。假定声源的初始方位已知,初始角速度为0,用本文第2节提出的算法和文献[4]中的算法做跟踪实验,画出水平角和仰角方向的声源轨迹,分别如图3和图4所示,其中虚线为真实声源轨迹,实线为估计的轨迹。由图3可见,本文提出的跟踪算法估计的方位与真实方位拟合得较好。图3的结果也验证了前文的分析,声源的水平角变化范围较大,而仰角的变化范围要小得多。用文献[4]中的算法估计出三维声源轨迹然后转换成方位估计,如图4所示。对照图3和图4可见,本文提出的跟踪算法比传统算法的方位估计精度要高。

表1给出了上述两种算法的水平角,仰角以及总的方位跟踪均方根误差。由该表可见本文提出的算法比传统算法的方位跟踪均方根误差有明显减少。

此外,由表1可以看出仰角方向的跟踪精度比水平角方向的要低,从图3和图4也可以得出同样的结论。这是因为实验中采用的是平面阵,该阵型对水平角有较高的分辨率,而对仰角则不敏感[12]。在麦克风阵列的实际应用中,如视频会议,声源的仰角变化很小,水平角变化很大,因而跟踪算法的主要作用是提供水平角的估计值。均匀圆阵在360o 的水平角范围内具有均匀的角度分辨率,非常适合这种应用场景。

4 结 语

传统的跟踪算法基于大型麦克风阵列对声源进行三维跟踪。而麦克风阵列的很多实际应用只需要估计声源的方位。本文提出一种基于粒子滤波的声源方位跟踪算法。该算法在Langevin方程的基础上建立移动声源的方位变化模型,采用SRP-PHAT作为定位函数,运用粒子滤波对移动声源进行方位跟踪。实验表明,本文提出的算法可在真实环境中实现随机走动说话人的方位跟踪,并得到了较高的方位跟踪精度。该算法适合于小型阵列、室内环境中的说话人方位跟踪。

参考文献

[1]Zhang Cha,Florenciol D,Ba D E,et al.Maximum likelihood soundsource localization and beamforming for directional microphone arraysin distributed meetings[J].IEEE Transactions on Multimedia,2008,10(3):538-548.

[2]Avokh A,Abutalebi H R.Speech enhancement using linearly con-strained adaptive constant directivity beam-formers[J].Applied A-coustics,2010,71(3):262-268.

[3]Markovic I,Petrovic I.Speaker localization and tracking with a micro-phone array on a mobile robot using von Mises distribution and particlefiltering[J].Robotics and Autonomous Systems,2010,58(11):1185-1196.

[4]Vermaak J,Blake A.Nonlinear filtering for speaker tracking in noisyand reverberant environments[C]//Proceedings of IEEE InternationalConference on Acoustics,Speech,and Signal processing(ICASSP),Salt Lake City,UT,USA,2001:3021-3024.

[5]Ward D B,Lehmann E A,Williamson R C.Particle filtering algo-rithms for tracking an acoustic source in a reverberant environment[J].IEEE Transactions on Speech and Audio Processing,2003,11(6):826-836.

[6]Liang Zhiwei,Ma Xudong,Dai Xianzhong.Robust tracking of movingsound source using scaled unscented particle filter[J].Applied A-coustics,2008,69(8):673-680.

[7]Talantzis F.An acoustic source localization and tracking framework u-sing particle filtering and information theory[J].IEEE Transactionson Audio,Speech,and Language Processing,2010,18(7):1806-1817.

[8]曹洁,李伟.基于正则化粒子滤波的说话人跟踪方法[J].兰州理工大学学报,2010,36(6):85-88.

[9]Dibiase J H.A high-accuracy,low-latency technique for talker locali-zation in reverberant environments using microphone arrays[D].Prov-idence:Brown University,2000.

[10]Arulampalam M S,Maskell S,Gordon N,et al.A Tutorial on particlefilters for online nonliner/non-Gaussian Bayesian tracking[J].IEEETransactions on Signal Processing,2002,50(2):174-188.

[11]Lathoud G,Odobez J M,Gatica-Perez D.AV 16.3:an audio-visualcorpus for speaker localization and tracking[R].Martigny:IDIAPResearch Institute,2004.