跟踪滤波算法仿真分析

2025-01-16

跟踪滤波算法仿真分析（精选7篇）

跟踪滤波算法仿真分析篇1

由于真实客观环境中噪声的特性常常是不平稳的, 因而自适应滤波在噪声对消中的应用有着广阔的发展前景。随着自适应噪声对消技术的发展, 应用于自适应噪声对消系统的算法也发生了很大的变迁。其中最常用的是LMS (最小均方算法) 和RLS (递归最小二乘算法) , 而LMS算法[1]以其算法简单, 运算量小, 实现容易等优点而得到了广泛应用, 但固定步长的LMS算法在收敛速率, 跟踪速率和权失调噪声要求方面存在矛盾;为了解决这一矛盾, 人们研究出了变步长的LMS算法。但基于变步长LMS的算法存在收敛特性和失调量受步长影响的缺点, 且最优步长不太容易确定, 因此NLMS, ELMS等关于LMS的改进的滤波算法得到了发展。RLS这种算法是对输入信号自相关的矩阵求逆, 进而不断递推估计实现权值更新的, 具有更快的收敛速度, 但这种算法计算起来较为复杂, 存储量大, 不适用于实时性要求高的场合;另外一个使用局限是, 输入信号自相关的矩阵求逆之后必须具有正定性, 否则会引起算法发散。

1 算法介绍

1.1 最小均方 (LMS) 算法

LMS算法是随机梯度算法家族中的一员, 简单性是它的一个显著特点, 而且它不需要计算有关的相关函数, 也不需要矩阵求逆运算, 因此它也是线性自适应滤波算法的参考标准[2]。

LMS算法采用的是一种瞬时估计, 即用n时刻的平方误差性能函数|e (n) |2作为瞬时均方误差ξ=E[e2 (n) ]的估值, 其实质是以当前输出误差、当前参考信号和当前权系数求得下个时刻的权系数。

LMS算法输出信号y (n) 、输出误差e (n) 及权系数W (n) 的计算公式为:

其中μ是控制自适应速度与稳定性的增益常数, 称为步长因子, 选择时, 应该综合考虑收敛速度和稳态误差的要求。

自适应滤波器收敛的条件是:

, λmax是输入信号的自相关矩阵的最大特征值。

LMS算法的优点是结构较为简单, 适应变化能力强, 但其则具有收敛的速度较慢的缺点。为了能适用于信号实时性处理的场合, 如何提高LMS这种算法的自适应速度就显得尤为重要。

局限LMS算法收敛这一要素的主要原因有:

1) 步长因子不能过大, 不然算法最终不收敛;

2) 收敛速度及均方误差不能兼得。

这两个原因都与步长有关, LMS算法中的步长是唯一能够控制算法迭代过程的参量, 必然是改进LMS算法性能的唯一着手点。

1.2 归一化最小均方 (NLMS) 算法

LMS算法是通过对梯度矢量各分量单个数据取样值的估计得到的, 没有进行平均, 才会使梯度估计中存在着噪声。NLMS引入变步长的迭代过程[3], 加快了收敛速度。

NLMS算法的输出信号y (n) 、输出误差e (n) 及权系数W (n) 的计算公式为:

此即归一化LMS算法, 其步长被输入信号的范数平方除, 因而较LMS算法具有更好的稳定性和收敛性。

1.3 递归最小二乘 (RLS) 算法

RLS算法是一个递归实现, 其收敛速率比一般的LMS滤波器快一个数量级, 因此它在线性自适应滤波器中应用非常广泛[4]。

RLS算法的输出信号y (n) 、输出误差e (n) 及权系数W (n) 的计算公式为:

其中, 增益矢量g (n) =C (n-1) X (n) /[λ+XT (n) C (n-1) X (n) ];C (n) 为自相关矩阵Rxx (n) 的逆矩阵, 其定义式为C (n) =λ-1[C (n-1) -g (n) XT (n) C (n-1) ], 且C (0) =δ-1I (I为单位矩阵, δ为小的正实数) ;常数λ是遗忘因子, 要求0<λ≤1。

RLS算法主要应用于系统辨识、自适应控制和自适应信号处理等领域。主要优点是收敛速度快, 因此在快速信道均衡、实时系统辨识和实际序列分析中得到广泛的应用[4], 其主要缺点是每次迭代计算量大。

2 自适应算法的MATLAB仿真

2.1 LMS算法的MATLAB仿真

输入为正弦信号与随机噪声的迭加, 随机噪声的幅值小于1, 在取不同步长情况下, LMS算法的误差函数曲线如图1所示, 误差函数图中纵坐标表示误差的大小。

2.2 NLMS算法的MATLAB仿真

1) 输入为正弦信号与随机噪声的迭加, 随机噪声的幅值小于1。NLMS均值曲线图仿真结果如图2所示。每个图中的横坐标都表示迭代次数, 学习曲线中纵坐标表示均方误差的大小。

2) 当输入信号为正弦函数与噪声的叠加时, LMS和NLMS的性能对比如图3所示。

2.3 RLS算法的MATLAB仿真

仿真采用多权的自适应横向滤波系统, 期望响应是一个经过滤波的高斯随机噪声, 采用RLS算法的自适应滤波学习曲线和矢量估计误差曲线如图4所示。

3 结果分析

1) LMS算法最大的优点是易于实现, 而且对有限寄存器长度造成的实现误差不敏感, 在实际生活和生产中应用较为广泛。

由图1三组不同情况下LMS算法的误差图的对比可以看出, 开始时误差比较类似于正弦函数, 随着自适应过程的进行, 误差越来越小并且随机性增大, 随着迭代次数的增加逐渐趋于零附近。由图 (b) 、图 (c) 比较可知, u越小, 收敛速度越慢, 但稳态误差较小;由图 (a) 、图 (b) 比较可知, 阶数k越小, 收敛越快, 但稳态误差较大。从而验证了在迭代收敛过程中, 误差函数随着迭代次数的增加逐渐趋于零, 学习曲线也趋于在附近小幅度波动, 甚至为零。但是LMS算法的收敛速度和其稳定性能是相互矛盾的;步长较大时收敛速度较快, 但其稳定性较差;步长较小时收敛速度较慢, 但其稳定性较好。

2) 由图2可知, u越大, 曲线收敛的越快, 越容易趋于零, 但曲线却更不光滑, 振荡较大, 符合NLMS算法的规律。另外, 在NLMS算法中, 当u太大时, 学习曲线反而会发散;克服了LMS算法的缺点, 算法本身可看成是一种变步长的自适应算法, 它的步长大小与输入信号的信噪比有关。

3) RLS算法在收敛速度和信号稳定性方面的性能都比LMS和NLMS算法良好, 收敛速度比LMS算法快一个数量级, 收敛性能与输入信号的频谱特性无关而且对信号的跟踪能力较强, 误差较小。但是RLS算法涉及到矩阵求逆, 计算复杂度很高, 所需的存储量极大, 不利于实时实现[5]。

4 结论

本论文主要介绍了三种常用的自适应算法:LMS、NLMS及RLS, 并通过MATLAB仿真, 从收敛性、误差函数和学习曲线等方面对这三种算法进行了简单的分析。结果表明, LMS算法易于广泛的应用, 但步长因子存在难以调和的矛盾;NLMS在收敛速度上有了明显的改进, 是对于LMS的优化;RLS收敛速度和稳定性都很好, 但计算量过大, 不利于大范围推广。

摘要：主要对自适应滤波算法展开了研究和讨论, 重点对LMS算法、NLMS算法以及RLS算法做了详细的说明和对比, 在算法原理、算法性能分析方面说明了各自算法的优越性。通过MATLAB仿真, 对每种算法的收敛性、学习曲线和误差分析等方面进行了分析。

关键词：自适应,噪声对消,LMS算法,NLMS算法,RLS算法

参考文献

[1]王海峰, 陈伟, 黄秋元.基于LMS算法自适应噪声抵消器的分析研究[J].计算机与数字工程, 2009, 37 (3) :85-87.

[2]秦彦平.基于DSP地下漏水定位系统的自适应滤波器设计[D].内蒙古大学, 2010.

[3]文静, 文玉梅, 李平.基于噪声白化准则的自适应噪声抵消方法[J].仪器仪表学报, 2010, 31 (8) :1693-1698.

[4]肖林.LMS和RLS算法的研究与实现[J].中山大学研究生期刊, 2010, 31 (2) :73-81.

[5]余慧.基于DSP反馈自适应噪声抵消器研究与设计[D].湖南大学, 2010.

跟踪滤波算法仿真分析篇2

关键词：交互式多模型,粒子滤波,声源跟踪

在机动目标跟踪领域,交互式多模型(interact-ing multiple model,IMM)[1,2]是在第一代多模型算法的基础上引入模型之间的交互而得出的。经典IMM算法采用卡尔曼作为滤波器,而卡尔曼滤波器只能处理线性系统、高斯噪声情况下的跟踪,这就导致IMM算法应用的局限性。因为粒子滤波算法 (particle filter,PF)[3]解决了非线性系统、非高斯噪声情况下的跟踪问题,所以该算法成为近年来IMM研究的热点之一。

粒子滤波是基于蒙特卡罗仿真的递推最优贝叶斯估计。其基本思想是:通过寻找一组在状态空间中传播的随机样本对概率密度函数进行近似,根据样本均值及其对应的权值计算状态变量的估计值[4,5]。

为了将IMM算法应用在非线性、非高斯的环境下,Boers和Driessen[6]将IMM与PF结合,构建了交互式多模型粒子滤波(interacting multiple model parti-cle filter,IMMPF)。IMMPF算法在模型的交互和融合方面与原有的IMM算法相同,改进的地方在于滤波器由原来的卡尔曼替换为粒子滤波。相比于IMM算法,IMMPF算法能够在非线性、非高斯的跟踪情况下得到状态输出的最优解。

1基本的IMMPF算法

系统的状态转移方程和观测方程可以表示为

式中,X(k) 表示k时刻模型mk的状态向量;F和G分别为状态转移矩阵和过程噪声输入矩阵;W表示过程噪声。对于不同的模型,状态转移矩阵F 、过程噪声输入矩阵G以及过程噪声W是不相同的。 Z(k) 表示相应的观测向量;H表示观测矩阵;V表示观测噪声。

假设系统中模型个数为M ,mk= 1,2,⋯M 。系统模型之间的转移概率由马尔科夫链表示,即

一般来说,在整个声源运动过程中pij是恒定不变的。若已知初始状态X( )0 ,初始模型概率以及各个时刻观测值Z(k) ,就可以估计出k时刻的状态估计X(k) 。

假设k - 1时刻各个模型的模型概率为第i个模型在k - 1时刻的状态估计和相应的协方差分别表示为Xi(k - 1) 和Pi(k - 1) 。交互式多模型粒子滤波算法主要分为以下四步。

(1)输入交互

k - 1时刻模型的交互概率为

式中,归一化因子为

交互后模型j在k - 1时刻的输入为

交互后模型j在k - 1时刻状态变量的协方差为

(2)模型匹配滤波

预测k时刻粒子状态可表示为

式中,Fj为模型j的状态转移矩阵;Gj和Wj分别为模型j的过程噪声输入矩阵及过程噪声。

预测观测值如下

Hj为模型j的观测矩阵。

粒子的残差为,且

其中,R是观测噪声的协方差。

归一化后的权值为

k时刻模型j的状态估计为

k时刻模型j协方差估计为

(3)模型概率修正

预测观测值的均值:

残差协方差为

(4)交互输出

2改进的IMMPF算法

为了提高IMMPF算法的跟踪精度、降低算法的运算量,相关学者从不同的角度对IMMPF进行了改进。改进算法大致可以分为三类:

(1)基于马尔科夫转移概率的修正

经典的IMMPF或IMM算法一般假设马尔科夫转移概率为固定值且是一阶的。一阶的马尔科夫转移概率容易造成滤波的精度不高。为了提高跟踪的精度,JIAN LAN等[7]将二阶马尔科夫转移概率引入到IMMPF算法中,文献通过实验证实采用二阶马尔科夫转移概率的IMMPF算法与原有算法相比,跟踪精度更高。封普文等[8]为了打破转移概率为固定值的限制,利用后验信息对马尔科夫转移概率矩阵进行修正,提出一种马尔科夫转移概率矩阵在线更新算法。IMMPF算法的模型集是有限的,在跟踪强机动目标时容易出现发散的现象。文献[9]提出一种机动性检测的PTHMM算法。该算法将目标的机动性作为观测量,马尔科夫转移概率矩阵作为隐藏状态,建立隐马尔科夫模型,利用维特比算法计算出模型转移概率权值,进而动态修正转移概率矩阵。

马尔科夫转移概率矩阵是交互式多模型粒子滤波很重要的一个参数,其精确与否将直接决定交互后每个滤波器输入状态变量的精确度,进而影响整个算法的精确度。

(2)基于模型概率的修正

模型概率作为决定最后输出结果的一个参量, 其对最终的输出结果也起着至关重要的作用。模型概率类似于粒子滤波中权值,对并行运算的每个滤波器的输出结果加权融合得到最终的输出结果。对模型概率进行修正也能提高算法的跟踪精度。

在IMMPF算法中,模型概率的计算只是利用了新息和模型概率的预测值,没有利用当前时刻的状态协方差。文献[10]综合利用状态协方差的信息、新息和模型概率的预测值等信息,提出基于模型概率修正的交互式多模型算法。朱军祥[11]基于模糊推理算法对模型概率进行了修正。修正后的算法将较小的模型概率设置成零,使与当前运动状态匹配程度更高的模型概率变大,进而降低模型之间的竞争, 提高跟踪的精度。

(3)滤波算法的改进

IMMPF算法在解决非线性、非高斯问题时,能够获得最优解。但当目标处于线性运动时,IMMPF算法的跟踪精度要低于IMM算法。文献[12]对滤波部分作了改进,其在滤波部分综合使用了卡尔曼滤波和粒子滤波。当状态方程是线性时滤波器匹配卡尔曼滤波器,当状态方程是非线性模型匹配粒子滤波。虽然粒子滤波在解决非线性、非高斯问题时具有无可比拟的优势。文献[13]为解决粒子滤波退化问题,提出一种改进的残差重抽样算法。算法不再单独求取每个粒子的复制次数,而是对粒子权重与总的粒子数的乘积累积求和后再取整,也就是求得累积粒子复制次数,最后输出粒子,保证重采样前后粒子数目不变。文献[14]通过改变目标状态方程的似然函数来调整粒子的权值,使得采样粒子更有效。文献[15]将无迹粒子滤波运用到IMM算法。无迹粒子滤波的重要性抽样函数考虑了当前时刻的量测信息,采样得到的粒子更接近真实后验概率密度函数,从而使滤波结果更加精确。

3IMMPF算法在声源跟踪的仿真结果与分析

声源跟踪实质上是一种非线性滤波,运动模型是影响跟踪效果很重要的一个因素。IMMPF的运动模型可以包含多种运动形式,因此可以跟踪随意运动的目标。声源跟踪作为目标跟踪的一种应用场景,将IMMPF算法应用在声源跟踪当中一方面丰富了声源跟踪的理论,另一方面也拓展了IMMPF算法的应用领域。下文的论述是对IMMPF算法在声源跟踪的应用进行了仿真实验。

3.1实验仿真环境

实验在一个5 m×4 m×3 m的房间内进行,房间的混响通过IMAGE模型仿真获得。实验中选取房间混响为300 ms,信噪比为30 d B的一段语音信号进行仿真分析。

实验采用匀速、左转弯和右转弯三种模型对声源的运动轨迹进行模拟仿真。实验中观测时间间隔T = 0.064 s,初始模型概率均为1 3 ,模型转移概率矩阵,粒子数N = 200 。匀速模型、左转弯以及右转弯三个模型对应的过程噪声协方差Q1= 0.01,Q2= 0.1,Q3= 0.1 ,观测模型噪声协方差R = 0.01 。三种模型初始状态协方差分别为、P1(0)=diag(10- 6× [400,100,400,100])、P2(0)=diag(10-6×[100,100,100,100])、P3(0)=diag(10- 6× [100,100,100,100])。初始状态设定为[0.8; 0.4; 0.2; 0.3] ,转弯速率为1.13°/s。

运动轨迹:

整个声源运动时间持续100 s,主要经历五个阶段

(1)0~20 s匀速直线运动;

(2)21~45 s声源以1.13°/s的速率作左转弯运动;

(3)46~60 s目标继续作匀速运动;

(4)61~85 s声源以1.13°/s的速率作右转弯运动;

()匀速运动。

3.2实验仿真结果分析

实验仿真结果如图1所示。由图1可以看出, IMMPF算法在匀速运动阶段能够得到较好的跟踪效果。当声源处于转弯运动时,跟踪结果与真实轨迹稍有偏差。转弯模型是影响跟踪精度的一个很重要的参量,对转弯模型的优化有助于提高跟踪的精度。

4结论

IMMPF算法融合了IMM和PF各自的优点,在解决非线性、非高斯情况下的状态估计具有明显的优势。与此同时也要看到,IMMPF算法作为最近几年刚被提出的算法,其算法本身还不够成熟,在声源跟踪的应用也只是初步的尝试[16]。如要完善IMMPF算法在声源跟踪的应用,仍有许多问题亟待解决。

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

基于麦克风阵列的声源定位与跟踪技术可广泛用于视频会议[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帧时的状态矢量可表示为:

$α_{t} = [θ_{t}, ϕ_{t}, \bar{θ}_{t}, \bar{ϕ}_{t}]^{Τ}$ (3)

其中[θt,ϕt]T表示声源方位,[ $\bar{θ}$ t, $\bar{ϕ}$ 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时刻的声源状态估计 $\hat{α}$ t。贝叶斯滤波问题的解包括预测和更新两个步骤。假设t-1时刻的后验概率密度函数p(αt-1|Y1:t-1)是已知的,则t时刻的后验概率密度函数p(αt|Y1:t)可通过以下两个方程得到:

$p (α_{t} | Y_{1 : t - 1}) = \int \begin{matrix} p (α_{t} | α_{t - 1}) \end{matrix} p (α_{t - 1} | Y_{1 : t - 1}) d α_{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{θ}_{t} = a_{θ} \dot{θ}_{t - 1} + b_{θ} F_{θ}$ (9)

$θ_{t} = θ_{t - 1} + Δ Τ \dot{θ}_{t}$ (10)

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

$b_{θ} = \bar{v}_{θ} \sqrt{1 - a_{θ}^{2}}$ (12)

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

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

$f (l ‚ X_{t}) = \sum_{l = 1}^{Μ} \sum_{m = l + 1}^{Μ} \hat{R}_{l m} [τ_{l m} (l)]$ (13)

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

$\hat{R}_{l m} (τ) = \frac{1}{Κ} \sum_{k = 0}^{Κ - 1} \frac{X_{l} (k) X_{m}^{*} (k)}{| X_{l} (k) X_{m}^{*} (k) |} e^{j ω τ}$ (14)

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

$τ_{l m} (l) = \frac{ζ^{Τ} (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{α}_{t}^{(i)} = g (α_{t - 1}^{(i)}, u_{t})$ ;② 通过定位函数将原始数据转变成目标状态观测值,Yt=f(l,Xt);

③ 形成似然函数, $p (Y_{t} | \tilde{α}_{t}) = F (Y_{t}, \tilde{α}_{t})$ ;

④ 根据似然函数计算新的粒子权重, $\tilde{w}_{t}^{(i)} = p (Y_{t} | \tilde{α}_{t}^{(i)})$ ,然后将权重归一化, $w_{t}^{(i)} = \tilde{w}_{t}^{(i)} \cdot (\sum_{i = 1}^{Ν} \tilde{w}_{t}^{(i)})^{- 1}$ ;

⑤ 计算当前时刻声源方位估计, $\hat{l}_{s} (t) = \sum_{i = 1}^{Ν} w_{t}^{(i)} l_{\tilde{α}}^{(i)}$ ;

⑥ 按照权重,从现有粒子组 ${\tilde{α}_{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)来评价,定义为:

$Ε_{R Μ S E} = \sqrt{\frac{1}{Κ_{s}} \sum_{t = 1}^{Κ_{s}} l_{s} (t) - \hat{l}_{s} (t)^{2}}$ (18)

其中,ls(t)表示第t帧时刻声源的真实方位, $\hat{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.

跟踪滤波算法仿真分析篇4

关键词：粒子滤波,多变量拉普拉斯分布,颜色直方图,重尾问题

粒子滤波作为一种后验概率的求解方法, 通过蒙特卡罗方法实现递推贝叶斯滤波, 在处理非高斯和非线性时变系统的参数估计和状态滤波问题方面具有独特的优势, 近年来已成为目标跟踪的一个强有力的工具[1-4]。但由于算法还处在发展阶段, 算法本身还不是很成熟, 仍存在着一些缺陷, 使得粒子滤波在实际应用中还存在着一些问题。在实际应用过程中, 发现在粒子滤波框架下采用颜色直方图作为目标特征, 当目标突然快速转向时, 跟踪效果往往不好, 甚至有时跟踪失败。视频运动目标的机动变化会使跟踪性能大大恶化, 尤其是当目标以不规则的速度进出某一场景, 产生复杂的交互运动和遮挡等问题的时候, 往往会导致目标跟踪失败。运动目标突然快速转向, 会产生重尾问题, 这个问题常常被很多研究者所忽视。多变量拉普拉斯分布可以很好地处理重尾问题。因此, 提出一种基于变换观测模型的粒子滤波跟踪算法, 提高视频目标的跟踪精度。

1 重尾问题

选取颜色直方图作为目标特征, 当视频目标在一段时间稳定运动, 突然发生偏转的时候, 就会出现重尾现象。图1给出了颜色直方图的2个重尾现象的实例。图1a是一个自然图像的对数直方图, 其中的黑色实线为跟踪目标的直方图结果, 灰色虚线为高斯分布的拟合结果, 从中可以看出黑色实线即目标的直方图中间部位是一个突起的尖峰, 而两边缓慢下降且拖尾较长。显然, 常用的高斯分布观测模型与目标的实际偏差较大, 必须寻求其他分布来更好地逼近目标的真实分布;图1b是一张拍摄图像经过小波变换后取其第二层分解的低频子带系数的直方图, 重尾现象更加明显。如何才能逼近目标快速突变转向时的真实分布, 为了解决这个问题, 在这里引入多变量拉普拉斯分布。

2 多变量拉普拉斯分布

一般的多变量拉普拉斯分布可以表示为[5]

式中, X, Y∈Rd, X～N (0, σ2Id) ;pZ (z) =exp (-z) 且z≥0, 记为Y～ML (μ, σ2) , 其概率密度函数为

式中, ;Km (x) 为改进的第二类贝赛尔函数;m为贝赛尔函数的阶数;d为变量y的维数。由贝赛尔函数渐近线方程, 当|x|→∞时, 可以得到

图2是一个多变量拉普拉斯分布图, 结合式 (3) 和图2可以看出, 与高斯分布相比, 多变量拉普拉斯分布能够更好地抓住重尾特性, 非常适合对重尾问题进行处理。

3 观测模型的建立及变换准则

假设目标的状态为X= (x, y, ) T, 目标的运动模型为

式中, 为系统噪声。假设候选目标的中心位于l= (x, y) 处, 目标的颜色信息可以用直方图来表示[6]。假设目标参考模板的颜色分布为则候选目标与目标参考模板的相似程度可以用Bhattacharyya系数来衡量:, 两者之间的差异可以用Bhattacharyya距离来衡量:, 则k时刻处于状态的一个粒子, 其颜色信息的观测似然函数为

其中, λ为归一化常数, 通过试验取0.01。式 (4) 是用于目标在一段时间内稳定运动情形的观测似然函数。当目标突然转向时, 将出现重尾现象, 结合式 (2) , 令μ为归一化常数, 重尾问题的颜色信息观测似然函数为

所谓变换观测模型, 就是采用2个观测模型根据跟踪的实际进行变换, 当目标稳定运动时, 采用式 (4) 的观测似然函数, 而当目标快速突变转向时, 采用式 (5) 的观测似然函数。接下来的关键问题就是如何判断何时采用那种观测模型, 给出观测模型的变换准则。令为k时刻的N个粒子, 每个粒子的权值分别为, 则k时刻目标估计值为。k时刻目标估计的方差定义为

其中diag (M) 表示由矩阵M对角线上的元素构成的向量。取k时刻N个粒子标准差的均值即

将σk, m的值与设定的阈值γ作比较, 观测模型变换的规则如图3所示。如果σk, m>γ, 则采用式 (4) 所示的观测似然函数, 否则, 采用式 (5) 所示的似然函数。在这里阈值γ设定为目标矩形短半轴的长度。

4 算法步骤

综上所述, 基于变换观测模型的粒子滤波跟踪算法的详细步骤如下:

(1) 初始化:选择目标跟踪区域, 获得该矩形区域的参考颜色直方图, 根据已知的先验信息p (x0) , 初始化粒子状态{x0i}Ni=1;

(2) 选择式 (4) 作为观测似然函数;

(3) 根据目标运动方程进行状态传播, 得到下一时刻的新粒子;

(4) 进行重要性采样, 计算颜色直方图分布, 根据直方图距离, 结合观测似然函数, 计算粒子的重要性权值并归一化处理;

(5) 进行粒子重采样处理;

(6) 计算状态均值, 输出被估计状态;

(7) 计算粒子标准差的均值σk, m的值, 并与阈值γ比较大小, 若σk, m>γ, 则选取式 (5) 作为观测似然函数, 直接转到步骤 (3) , 否则返回步骤 (2) 。

5 仿真比较与分析

为了验证算法的跟踪性能, 分别对三组视频进行跟踪仿真。为了便于比较和说明问题, 分别采用基于粒子滤波的颜色直方图特征的跟踪算法和文中提出的算法进行跟踪测试。取400个粒子, 在PentiumⅣ2.8 GHz, 内存为1 GB的计算机上进行试验。视频的参数如表1所示, 各视频的跟踪仿真结果分别如图4、图5和图6所示。

由仿真结果可以看出:对于视频1, 粒子滤波的颜色直方图特征的跟踪算法在目标运动方向第一次变化时, 跟踪就有所偏差, 在目标运动方向第二次变化时偏差更大, 最终导致丢失目标;而这里提出的基于变换观测模型的粒子滤波跟踪算法可以较好地跟踪到目标。对于视频2, 粒子滤波的颜色直方图特征的跟踪算法在几次目标运动突然偏转时, 跟踪偏差较大, 在第51帧的时候跟踪发生错误, 跟踪到其他相近的目标, 而基于变换观测模型的粒子滤波跟踪算法在目标运动多次突然偏转的情况下, 一直跟踪效果较好。对于视频3, 目标在第14、206、321帧都突然转向, 基于粒子滤波的颜色直方图特征的跟踪算法在目标运动方向几次发生变化时, 都会丢失目标;而文中提出的基于变换观测模型的粒子滤波跟踪算法可以较好地跟踪到目标, 进一步验证了算法的有效性。

6 结束语

跟踪滤波算法仿真分析篇5

本文基于文献[1]提供的关于联合估计目标状态和系统误差的快速最大后验估计算法[3] (简称FMAP算法) , 采用数值仿真的方法, 研究算法的稳健性, 即最后通过变化算法输入参数, 以二维匀速运动目标 (两雷达同步量测) 为例, 采用数值仿真的方法分析算法对输入参数的依赖程度即稳健性。

1 数值仿真参数的设定

FM AP算法的输入参数主要包括三类:测距的标准差、测角的标准差、目标空间坐标白噪声的标准差。为对算法关于这三类输入参数的稳健性进行简单分析, 考虑固定其他参数及目标与量测性能参数, 变化单个输入参数若干次, 在同一坐标系中绘制每次的目标位置估计的RM SE曲线, 分析不同曲线的差异, 若变化趋势差异不大, 则从一定角度说明了算法的稳健性。这里以二维匀速目标二雷达情形为例进行仿真分析。目标运动参数与雷达参数的基本设置如下:

2 算法关于输入参数稳健性的仿真结果

2.1 算法对测距标准差的稳健性

考虑分别将雷达1的距离标准差依次调整为10, 20, …, 80, 固定其余设置, 仿真计算得到8条RMSE曲线的图形如下:

2.2 算法对测方位角标准差的稳健性

考虑分别将雷达1的方位角标准差依次调整0.05, 0.10, …, 0.45, 固定其余设置, 仿真计算得到9条RM SE曲线的图形如下:

2.3 算法对目标空间坐标标准差的稳健性

考虑分别将雷达1的距离标准差依次调整为10-11, 10-10, …, 103, 固定其余设置, 仿真计算得到8条RMSE曲线的图形如下:

仿真结果表明, 随着测量次数的变化, 输入参数的适度变化不会造成RMSE曲线变化趋势的变化, 甚至在曲线尾部取值偏离不会太远, 这说明若控制算法的输入参数在一定范围内取值, 目标位置的估计不至于出现较大的异动。

参考文献

[1]Zhou Lin, Pan Quan, Liang Yan, et al.Joint Estimation of State and Sensor Systematic Error in Hybrid System[J].2012.15th International Conference on Information Fusion, 2012.

[2]Helmick R E, Conte J E, Hoffman S A, et al.One-step fixed-lag interacting multiple model (IMM) smoothing for alignment of asynchronous sensors[J].in Proc SPIE Signal and Data Processing of Small Targets, 1994.

跟踪滤波算法仿真分析篇6

但是粒子滤波在迭代过程中会发生粒子退化,Doucet等人提出选用“好的”重要概率密度函数来解决粒子退化问题[5]。常用的方法是采用扩展卡尔曼滤波(extended kalman filter,EKF)、无迹卡尔曼滤波(unscented kalman filter,UKF)及容积卡尔曼滤波(cubature kalman filter,CKF)来重新设计重要性密度函数,由此得到扩展粒子滤波(extended particle filter,EPF)、无迹粒子滤波(unscented particle filter,UPF)及容积粒子滤波(cubature particle filter,CPF)算法。文献[6]中结合EKF提出EPF算法,但是EKF算法本身在处理非线性问题时,是依靠泰勒展开来近似线性化,精度上的偏差较大,容易发散,且计算复杂,在强非线性系统中应用较少[7]。文献[8]中提出了在闪烁噪声下使用UPF算法对系统进行滤波,其利用的UKF算法需要选取2n+1个sigma点计算系统的均值及方差,对参数设计要求较高,且在系统状态超过三维时容易出现发散[9,10,11]。由Arasaratnam等提出的CKF算法[12],它利用2n个等权值容积点代替复杂积分运算,且不用调整滤波参数,较UKF有着更高的精度与稳定性。

本文在闪烁噪声环境下,利用CKF设计粒子滤波的重要性密度函数[13],给出了CPF算法的简要过程,仿真验证了CPF算法的跟踪性能,并与PF和UPF算法作比较。

1 闪烁噪声

闪烁噪声是目标跟踪中一种典型的非高斯量测噪声。雷达跟踪目标时,由于复杂目标不同部位相对于雷达散射强度和相位的随机变化,引起回波波前面畸变产生的测量误差被称为闪烁噪声。闪烁噪声最显著的特点是长拖尾[14],中心部分则近似于高斯分布,如图1所示。远而小的目标发生的闪烁通常可以不考虑,但当跟踪大目标且距离又比较近时,闪烁噪声对跟踪效果有很大影响。

描述拖尾分布的数学分布通常有t分布、均匀分布、大方差高斯分布和拉普拉斯分布等。而描述闪烁噪声用的较多的方法是高斯分布和拉氏分布的合成,其密度函数表达式为

式(1)中ε∈[0,1]为闪烁概率,fg(x),fl(x)分布分别表示高斯分布和拉氏分布的密度函数。

2 标准粒子滤波算法

非线性离散系统表示如下:

式(2)中xk为k时刻状态向量,zk为量测向量,f(·)为状态函数,h(·)为量测函数,Qk-1和Rk分别为独立分布的系统噪声和量测噪声方差矩阵。

标准粒子滤波的算法流程如下。

Step1:k=0,由先验概率p(x0)进行随机采样,得到等权值粒子集{xi0}Ni=1。

Step2:在k时刻,更新粒子权值:

归一化:

状态变量和协方差估值为

Step3:重采样得到新的粒子集合。

Step4:经过状态方程f传播得到新的粒子xik+1/k。

Step5:另k=k+1,转到Step2。

传统的粒子滤波在滤波过程中没有使用含有最新量测信息的数据,产生了很大的权重方差,不能较好地逼近后验概率[15]。

3 容积粒子滤波算法

利用CKF算法重新设计粒子滤波的重要性密度函数,即得到CPF算法,其具体流程如下。

Step1:设粒子数为N,通过先验分布p(x0)采集粒子样本,得到初始粒子集{x0i,i=1,2,…,N}。初始粒子集的状态和协方差为

Step2:使用CKF更新每个粒子,设系统状态维数为n。

时间更新:

量测更新:

Step3:采样:

更新粒子权值并归一化:

Step4:采用MCMC方法重采样。获得一组等权值粒子集。

Step5:计算状态变量和协方差矩阵的估值。

4 仿真分析

设雷达处于坐标原点,目标在平面内做匀速转弯(CT)运动,状态方程为

式(26)中状态变量xk=[xk,vxk,yk,vyk],Φ为状态转移矩阵,Γ为状态噪声驱动阵,Wk-1为状态噪声,其方差为Q=diag(2,2)。设目标观测方程为

式(27)中

式(28)中rk为斜距,θk为方位角,Vk为观测噪声,斜距的观测噪声为闪烁噪声,闪烁概率ε取0.05。方位角观测噪声是高斯噪声,其均值为零,标准差为0.02 rad。采样粒子数为100,采样间隔T=1 s,仿真时间为100 s。给定初始状态x0=[500,100,1 000,5]T,初始协方差矩阵P0|0=diag(10,0.1,10,0.1),用PF,UPF与CPF进行对比仿真。

匀速转弯模型的状态转移矩阵为

。进行100次蒙特卡洛仿真,仿真结果如图2、图3、表1。

具体运算结果如表1所示。

(1)从图2、图3及表1看出,在闪烁噪声干扰下,PF算法的跟踪效果随时间推移逐渐变差,甚至可能出现发散情况;而UPF和CPF算法则保持了良好的跟踪性能,原因是UPF和CPF重新设计了重要性密度函数,包含了最新的量测信息,因此跟踪效果较好。

(2)从运行时间上看,CPF算法实时性比UPF好,用时仅是UPF算法的1/2左右。这是由于在进行非线性函数传播时,UKF算法需要选择2n+1个sigma点,CKF算法只需要2n个容积点。

(3)根据表中的数据可以看出CPF的跟踪误差约等于PF和UPF的1/5和1/2,这是因为在设计重要性密度函数时,UKF算法对参数设定要求高,选择不当会引起自协方差非正定,致使滤波性能下降。

5 结束语

本文针对非线性非高斯闪烁噪声下的目标跟踪问题,首先建立了闪烁噪声的数学模型;其次介绍了传统的粒子滤波,并指出其不足之处;然后利用CKF算法重新设计了重要性密度函数得到CPF算法;最后在闪烁噪声下检验CPF算法性能,并与PF及UPF算法作对比,仿真结果表明,在闪烁噪声下,CPF算法较PF和UPF算法有着更好的跟踪精度与稳定性,且相比于UPF算法实时性有所提高。

摘要：针对闪烁噪声下非线性非高斯系统的目标跟踪问题,首先建立了闪烁噪声的数学模型;然后分析了传统粒子滤波算法的优劣点,在此基础上,引入容积卡尔曼滤波算法,重新设计粒子滤波的重要性密度函数,提出用容积粒子滤波算法来跟踪目标。最后进行了仿真分析与对比。仿真结果表明,闪烁噪声条件下,容积粒子滤波算法的跟踪误差分别是传统粒子滤波算法和无迹粒子滤波算法的1/5和1/2,有更高的跟踪精度;而运行时间仅是无迹粒子滤波算法的1/2,且跟踪稳定性更好。

跟踪滤波算法仿真分析篇7

滤波和预测是目标跟踪的重要组成部分, 其目的是估计当前和未来时刻目标的运动状态, 包括位置、速度、加速度等。跟踪滤波器的种类很多, 其中最常用的是卡尔曼滤波器[1,2,3,4,5]。由卡尔曼滤波器的工作过程可知, 某个时刻目标状态的更新值等于该时刻的预测值再加上一个与增益有关的修正项。而要计算增益, 就必须计算协方差的一步预测、新息协方差和更新协方差。因而在卡尔曼滤波中增益的计算占了大部分的工作量。为了减少计算量, 就必须改变增益矩阵的计算方法。为此人们提出了常增益滤波器, 此时增益不再与协方差有关, 因而在滤波过程中可以离线计算, 这样就大大减少了计算量。

已有的常增益滤波器如α-β、α-β-γ的目标状态向量中只包含位置、速度或加速度分量。随着目标机动性的增强, 这类滤波器的跟踪性能有所下降。针对这些情况, 有人提出了改进的α-β滤波算法[6,7,8]和α-β-γ滤波算法[9,10], 使它能够跟踪含有加速度噪声的机动目标, 但它的状态向量也只包含速度或加速度分量, 且适用范围有限。为了更加准确地估计目标的运动状态, Kishore等人提出了基于卡尔曼滤波的包含加加速度状态分量的跟踪模型, 即Jerk模型[11]。此模型对高机动目标的跟踪性能有明显提高, 但计算量过大, 耗时较长。基于此, 本文给出了一种针对匀加加速度运动目标的包含位置、速度、加速度和加加速度分量的常增益滤波算法——α-β-γ-δ滤波模型。

通过对α-β、α-β-γ、α-β-γ-δ滤波模型的理论和实验分析可知, 它们分别对匀速、匀加速和匀加加速度运动目标有很好的跟踪精度, 但现实中完全做上述三种运动方式的目标并不多, 而大部分目标都可看作是上述三种运动方式的组合。因此本文针对单一滤波器对目标机动性发生改变时, 滤波精度降低, 甚至出现目标丢失的现象, 提出了将α-β、α-β-γ、α-β-γ-δ滤波器组合起来使用的思想。在滤波参数不变的情况下, 可以得到更高的跟踪精度, 扩大了常增益滤波器的运用范围。

1 α-β-γ-δ模型

在时常系统分段常数白色噪声过程[1]中 (即在每一个采样周期内加加速度是常数, 并且在各周期之间是独立的) , 可设目标的状态方程为:

X (k+1) =FX (k) +Γv (k) (1)

其中, v (k) 是标量零均值白噪声序列, 协方差为:

E[v (k) v (j) ]=σ $_{v}^{2}$ δkj (2)

令关于坐标x的状态向量为:

$X = [x \dot{x} \ddot{x} \overset{\dots}{x}]^{´} (3)$

其中X、 $\dot{x}$ 、 $\ddot{x}$ 和 $\overset{\dots}{x}$ 分别表示k时刻目标位置、速度、加速度和加加速度;状态转移矩阵、过程噪声分布矩阵分别为:

$F = (\begin{array}{l} 1 Τ Τ^{2} / 2 Τ^{3} / 6 \\ 0 1 Τ Τ^{2} / 2 \\ 0 0 1 Τ \\ 0 0 0 1 \end{array}) (4)$

Γ=[T3/6 T2/2 T 1 ]′ (5)

T为采样周期。

过程噪声的协方差为:

$Q = Γ σ_{v}^{2} Γ^{'} = (\begin{matrix} Τ^{6} / 36 & Τ^{5} / 12 & Τ^{4} / 6 & Τ^{3} / 6 \\ Τ^{5} / 12 & Τ^{4} / 4 & Τ^{3} / 2 & Τ^{2} / 2 \\ Τ^{4} / 6 & Τ^{3} / 2 & Τ^{2} & Τ \\ Τ^{3} / 6 & Τ^{2} / 2 & Τ & 1 \end{matrix}) σ_{v}^{2} (6)$

假定只有位置量测信息可以应用, 则量测方程为:

Z (k) =HX (k) +w (k) (7)

其中

$Η = [\begin{matrix} 1 & 0 & 0 & 0 \end{matrix}] (8)$

量测噪声w (k) 为标量零均值白噪声序列, 其方差为:

$E [w (k) w^{'} (j)] = R δ_{k j} \overset{Δ}{=} σ_{w}^{2} δ_{k j} (9)$

对于时常系统的状态估计协方差P (k|k) , 在一定的条件下将收敛到一个稳态值[1]。这样可设状态估计协方差矩阵分量的稳定值为:

$\lim_{k \to \infty} Ρ (k | k) = \lim_{k \to \infty} Ρ (k + 1 | k + 1) = [p_{i j}] (10)$

一步预测协方差分量的稳定值为:

$\lim_{k \to \infty} Ρ (k + 1 | k) = [m_{i j}] (11)$

其中[pij]、[mij]均为对称阵。

而对于滤波增益, 其符号记为:

$\begin{array}{l} Κ = \lim_{k \to \infty} Κ (k + 1) \overset{Δ}{=} [\begin{matrix} g_{1} & g_{2} & g_{3} & g_{4} \end{matrix}]^{'} \\ \overset{Δ}{=} [\begin{matrix} α & β / Τ & γ / Τ^{2} & δ / Τ^{3} \end{matrix}]^{'} (12) \end{array}$

新息协方差的表达式为:

$S = Η (\begin{matrix} m_{11} & m_{12} & m_{13} & m_{14} \\ m_{21} & m_{22} & m_{23} & m_{24} \\ m_{31} & m_{32} & m_{33} & m_{34} \\ m_{41} & m_{42} & m_{43} & m_{44} \end{matrix}) Η^{'} + R = m_{11} + σ_{w}^{2} (13)$

其中, 使用符号R=σ $_{w}^{2}$ 。滤波增益可按下式计算:

$\begin{array}{l} Κ = (\begin{matrix} m_{11} & m_{12} & m_{13} & m_{14} \\ m_{21} & m_{22} & m_{23} & m_{24} \\ m_{31} & m_{32} & m_{33} & m_{34} \\ m_{41} & m_{42} & m_{43} & m_{44} \end{matrix}) Η^{'} S^{- 1} \\ = (\begin{matrix} \frac{m_{11}}{m_{11} + σ_{w}^{2}} & \frac{m_{12}}{m_{11} + σ_{w}^{2}} & \frac{m_{13}}{m_{11} + σ_{w}^{2}} & \frac{m_{14}}{m_{11} + σ_{w}^{2}} \end{matrix})^{´} (14) \end{array}$

由式 (12) 、式 (14) 可得:

$g_{1} = \frac{m_{11}}{m_{11} + σ_{w}^{2}} (15)$

$g_{2} = \frac{m_{12}}{m_{11} + σ_{w}^{2}} (16)$

$g_{3} = \frac{m_{13}}{m_{11} + σ_{w}^{2}} (17)$

$g_{4} = \frac{m_{14}}{m_{11} + σ_{w}^{2}} (18)$

协方差矩阵的更新方程变为:

$\begin{array}{l} [p_{i j}] = (Ι - Κ Η) [m_{i j}] \\ = (\begin{matrix} (1 - g_{1}) m_{11} & (1 - g_{1}) m_{12} & (1 - g_{1}) m_{13} & (1 - g_{1}) m_{14} \\ m_{21} - g_{2} m_{11} & m_{22} - g_{2} m_{12} & m_{23} - g_{2} m_{13} & m_{24} - g_{2} m_{14} \\ m_{31} - g_{3} m_{11} & m_{32} - g_{3} m_{12} & m_{33} - g_{3} m_{13} & m_{34} - g_{3} m_{14} \\ m_{41} - g_{4} m_{11} & m_{42} - g_{4} m_{12} & m_{43} - g_{4} m_{13} & m_{44} - g_{4} m_{14} \end{matrix}) (19) \end{array}$

其中I为单位方阵。

协方差矩阵的预测方程可重写如下:

P (k|k) =F-1[P (k+1|k) -Q] (F-1) ′ (20)

其中, 从 (4) 可得:

$F^{- 1} = (\begin{matrix} 1 & - Τ & Τ^{2} / 2 & - Τ^{3} / 6 \\ 0 & 1 & - Τ & Τ^{2} / 2 \\ 0 & 0 & 1 & - Τ \\ 0 & 0 & 0 & 1 \end{matrix}) (21)$

由式 (19) 和式 (20) 对应分量相等, 消去某些项后可解得关于α、β、γ、δ的四个方程:

$\frac{δ^{2}}{1 - α} = Τ Ι^{2} (22)$

γ2=2βδ (23)

γ2+24αγ+12αδ-12β2-24δ=0 (24)

γ2+96β+4αδ-48α2-48αβ-8δ-12β2=0 (25)

其中, $Τ Ι = \frac{Τ^{3} σ_{v}}{σ_{w}}$ 称为跟踪指标, α、β、γ、δ可数值地求解。

2 组合滤波器的设计

由于滤波器的残差反映了滤波器的滤波精度, 因此可以根据残差的大小来调整各种滤波器在组合滤波输出中的比重。具体做法如下:

设 $\hat{X}_{i} (k) ‚ d_{i} (k) ‚ i$ =1, 2, 3分别为α-β、α-β-γ、α-β-γ-δ滤波器的滤波输出和残差:

$d_{i} (k) = z (k) - \hat{X}_{i} (k | k - 1) i = 1 ‚ 2 ‚ 3 (26)$

其中 $\hat{X}_{i} (k | k - 1) ‚ i$ =1, 2, 3分别为α-β、α-β-γ、α-β-γ-δ滤波器的预测值。

记 ρ (k) =[1/d $_{1}^{2}$ (k) +1/d $_{2}^{2}$ (k) +1/d $_{3}^{2}$ (k) ]-1 (27)

则组合滤波输出 $\hat{X} (k)$ 可如下计算:

$\hat{X} (k) = \frac{ρ (k) \hat{X}_{1} (k)}{d_{1}^{2} (k)} + \frac{ρ (k) \hat{X}_{2} (k)}{d_{2}^{2} (k)} + \frac{ρ (k) \hat{X}_{3} (k)}{d_{3}^{2} (k)} (28)$

式中 $\frac{ρ (k)}{d_{1}^{2} (k)}$ 是第i个滤波器在组合滤波输出中所占的权值。若第i个滤波器的模型准确, 残差di (k) 就小, 那么它在组合输出中所占的比重就大。

3 仿真结果及分析

为了评价组合滤波器的性能, 对一种典型的目标运动形式进行了Monte Carlo仿真。仿真实验中, 采样周期T为0.5s, 目标机动时间为100s, 量测噪声的标准偏差σw为100m。目标在0～20s做匀速运动, 初始位置为1000m, 初始速度为200m/s;21～40s做匀加加速度运动, 加加速度为5 m/s3;41～60s做匀加速度运动, 加速度为40s时的值;61～80s做匀加加速度运动, 加加速度为-5 m/s3;81～100s做匀速运动, 速度为80s时的值。仿真的结果如图1和图2所示。

图1、图2是α-β、α-β-γ、α-β-γ-δ、组合滤波分别对目标进行跟踪的位置、速度RMS误差的比较。从图中可以看出α-β滤波对匀速运动目标有较高的跟踪精度, α-β-γ滤波对匀加速运动目标有比较高的跟踪精度, α-β-γ-δ滤波对匀加加速运动目标有比较高的跟踪精度。当目标运动方式发生改变时, 三种滤波器都会出现滤波精度下降, 而组合滤波能够比较好地克服这个问题, 能保持一个比较稳定的滤波精度, 且精度比较高。

4 结论

本文对机动目标跟踪的滤波算法进行了研究, 推导了一种四维常增益滤波器, 并把它与另两种常增益滤波器进行了组合, 得到了一种组合滤波器。通过仿真看出该组合滤波算法的滤波精度比较高、适用范围广, 且算法简单、实现容易, 特别适合工程上的应用。

参考文献

[1]Bar-Shalom Y, Fortmann T E.Tracking and Data Association[M].NewYork:Acadmic Press, 1988.

[2]何友, 修建娟, 张晶炜, 等.雷达数据处理[M].北京:电子工业出版社, 2006.

[3]周宏仁, 敬忠良, 王陪德.机动目标跟踪[M].北京:国防工业出版社, 1991.

[4]杨万海.多传感器数据融合及其应用[M].西安:西安电子科技大学出版社, 2004.

[5]康耀红.数据融合理论与应用[M].第2版, 西安:西安电子科技大学出版社, 2006.

[6]王红亮, 李枫, 赵亦工.一种机动目标跟踪的自适应α-β滤波算法[J].雷达科学与技术, 2007, 5 (4) :278-282.

[7]王晓伟, 金宏斌, 杨龙坡, 等.基于α-β和δ-ε滤波器的机动目标跟踪[J].现代雷达, 2006, 28 (7) :56-58.

[8]金宏斌, 徐毓, 钱琼芬.机动目标跟踪的改进α-β算法[J].空军雷达学院学报, 2003, 17 (2) :27-29.

[9]Sudano J J.Theα-β-ΓTracking filter with a noisy jerk as the maneuvermodel[J].IEEE Transactions on Aerospace and Electronics Systems, 1993, 29 (4) :578-580.

[10]Arcasoy C C.Theα-β-ΓTracking filter with a noisy jerk as the ma-neuver model:Frequency Domain Solution[J].IEEE Transactions onAerospace and Electronics Systems, 1996, 32 (3) :1170-1174.

【跟踪滤波算法仿真分析】推荐阅读：

组合滤波算法12-30