0 引言

粒子滤波是一种基于Monte Carlo思想的滤波算法.该算法将复杂积分转化为加权求和运算, 可有效估计非线性动态模型的隐状态, 目前已广泛应用于动态目标跟踪等领域.粒子滤波算法主要由重要性采样和重采样两部分组成, 其中重要性采样产生加权粒子集, 随着迭代步数增加会出现粒子退化现象.为了减弱该现象, Gordon等^[1]提出了重采样算法.该算法通过复制大权值粒子, 删除小权值粒子, 能够改善粒子退化问题, 但造成了粒子多样性匮乏^[2], 即产生粒子贫化问题.

为了解决粒子贫化问题, 针对重采样的改进算法一直是研究的热点^[3].最新的研究成果有:蒋蔚^[4]提出了一种基于SVM重采样的似然粒子滤波算法, 利用SVM重构后验概率密度来获得具有高多样性的粒子集, 改善粒子贫化问题; 张毅等^[5]提出了基于高斯分布重采样的RBPF-SLAM算法, 利用高斯分布散化高权重粒子, 得到新粒子集; Shen^[6]提出了一种智能粒子滤波算法IPF, 通过引入遗传算法中的交叉变异操作, 改善了粒子贫化问题.以上算法均需要预设完整的模型参数, 自适应性较差, 在目标系统出现状态突变时, 算法性能下降严重.

本文针对粒子滤波在跟踪非线性状态突变系统的隐状态时, 因粒子贫化出现的估计精度下降问题, 提出一种基于Student's t分布的自适应智能重采样粒子滤波算法.该算法通过引入Student's t分布的CDF作为参数转移方程, 自适应地导向生成大权值和小权值两个子集, 再对大小权值子集作自适应加权交叉和突变操作, 有效提升粒子的多样性.当系统出现状态突变时, 可自适应地保持粒子多样性, 从而达到提升算法精度和鲁棒性的目的.仿真实验的结果验证了本文所提算法的有效性.

1 粒子滤波算法分析

在非线性动态模型隐状态估计问题中, 首先针对目标问题构建动态模型, 该模型由如下状态转移方程和状态观测方程组成:

(1)

(2)

其中: f(x_k-1)和h(x_k)为非线性方程, N_k-1为状态转移噪声, R_k为状态观测噪声.

粒子滤波算法基于Monte Carlo^[7]学派的思想, 通过带权值的随机粒子集合{x_0:kⁱ, w_1:kⁱ}_N, 近似表示基于观测的后验概率密度函数p(x_0:k|z_1:k), 如下所示:

(3)

其中: k为迭代步数, i为粒子号, N为粒子总数.

然后, 通过后验概率密度函数估计目标状态的期望, 如下所示:

(4)

由于观测方程和状态方程的非线性, 无法对p(x_0:k|z_1:k)直接采样, 因此引入便于采样的重要性分布q(x_0:k), 如下所示:

(5)

其中可看作对重要性分布函数q(x_0:k)的加权, 权值用w_0:k表示为

(6)

因为在状态转移过程中, 第k步的目标状态只与第k-1步相关, 所以可将目标隐状态{x_0:k}看作Markov Chain.因此, 对式(6)作推导, 可得

(7)

(8)

将式(7)代入(8), 得到权值递推公式为

(9)

其中p(z_k|x_k)和p(x_k|x_k-1)均由动态模型得到, 并取q(x_k|x_k-1)=p(x_k|x_k-1).得到加权粒子集{x_kⁱ, w_kⁱ}_N的更新过程, 如下所示:

(10)

(11)

最后, 将权值w_kⁱ归一化后, 代入式(3)和(4), 得到后验概率密度的近似分布和目标隐状态的估计值.

在粒子滤波迭代中, 少数粒子权值变大, 趋向于1, 大量粒子的权值变小, 趋向于0^[8], 粒子集有效性降低.针对该问题, Gordon提出重采样算法.

重采样算法将粒子的权值统一为1/N, 使少数粒子权值过大的问题得到解决.但为了保存粒子集的特性, 重采样算法复制权值较大的粒子, 删除权值较小的粒子, 导致相同粒子增加, 产生粒子贫化^[9]现象.由于粒子滤波的核心思想为Monte Carlo加权求和, 粒子贫化直接导致算法精度下降.目前的新型重采样算法需预设大量参数, 对突变状态的自适应能力较差.针对该问题, 本文提出一种可根据目标状态动态自适应调整粒子集多样性的粒子滤波算法.

2 自适应智能重采样算法

当粒子权值和粒子状态通过式(10)和(11)更新后, 将粒子集分为包含较大权值的子集C_L和包含较小权值的子集C_S, 如下所示:

(12)

自适应分界值W_T的设计思路如下:

首先, 计算粒子有效采样尺度N_eff. N_eff可动态反映粒子集的有效性, 表示为

(13)

当N_eff偏大时, 粒子集有效性高.令C_L集合变大, 尽量保存粒子集原始状态, 使粒子集保持有效性.同时, 保留一部分小权值粒子, 以保持粒子多样性.当N_eff偏小时, 粒子集有效性低.令C_S集合变大, 使粒子集有效性较低时, 能够实现对更多小权值粒子的调整.因此, 选择Student's t分布的CDF函数作为参数转移方程, CDF特性如图 1所示.

由图 1可知, v=1时最符合要求, 因此将N_eff代入式(14), 得到自适应转移参数N_z为

(14)

然后, 将粒子的权值按大小递减顺序储存于中, 如下所示:

(15)

最后, 令第N_z号粒子的权值为分界值, 表示为

(16)

对小权值集合中的粒子进行交叉加权操作, 将权值大的粒子和权值小的粒子按照下式加权合并:

(17)

其中: x_kL^l ∈ C_L和x_kS^j ∈ C_S, x_kS^*j表示新生粒子.该操作可使小权值粒子向大权值粒子移动, 向高似然区域集中, 在重采样复制删除操作中, 保留更多的粒子, 从而达到提升粒子多样性的目的.另外, 由于更多的粒子处于高似然区域, 也使估计精度提高.权值α决定在交叉过程中大小权值粒子在新生粒子中所占的比例, 因此这里引入自适应算子, 依据粒子集的状态, 动态调整交叉权值α.

利用Student's t的CDF函数来获得自适应交叉权值α, 如下所示:

(18)

为了达到在加权交叉操作中, 粒子集有效性较小时, 不过度向大权值粒子集中, 有效性大时, 不过度拒绝加权交叉操作的目的, 结合图 1, 令参数v=0.1.

由Monte Carlo思想, 提升粒子集的多样性可以提升算法精度, 因此引入突变策略, 如下所示:

(19)

其中r_j为服从0到1均匀分布的随机变量.由于突变具有随机性, 为了防止过度改变粒子集, 将p_M预设为一个较小的常数0.3.

最后, 对调整后的新生粒子集进行重采样.由于自适应地调整了粒子集中的粒子分布, 保存了更多有效的粒子, 使粒子集多样性得到了提升.

3 算法流程

本文提出的基于Student's t分布的自适应重采样粒子滤波算法的完整步骤如下.

Step 1:初始化粒子集{x₀ⁱ}_N.

Step 2:对第k步进行隐状态估计, 利用式(10)和(11)更新粒子集, 得到加权粒子集合{x_kⁱ, w_kⁱ}_N.

Step 3:利用式(3)和(4)加权求和, 得到目标隐状态的估计值.

Step 4:利用式(13)得到粒子集的有效采样尺度N_eff.

Step 5:利用式(14) ~ (16)自适应地得到大小权值子集分界值W_T.

Step 6:将W_T代入式(12), 生成大权值粒子集C_L和小权值粒子集C_S.

Step 7:利用式(17)和(18)执行自适应加权交叉操作.

Step 8:利用式(19)执行突变操作, 最终得到新生粒子集合{x_k^*i}_N.

Step 9:将新生粒子集合{x_k^*i}_N代入式(11), 更新权值, 得到加权新生粒子集{x_k^*i, w_k^*i}_N, 执行重采样步骤.

Step 10:将{x_k^*i, w_k^*i}_N作为第k+1步的初始粒子集, 循环执行Step 2 ~ Step 10.

4 仿真实验

为了验证本文算法对非线性动态系统隐状态的跟踪性能和对非线性突变状态系统的跟踪性能, 采用一种具有强非线性的动态模型^[10], 其状态转移方程和状态观测方程如下所示:

状态转移模型为

(20)

状态观测模型为

(21)

其中:预设转移噪声为N_k~N(0, 4), 观测噪声为R_k~N(0, 4).

为了有效说明本文算法的特性, 选取SIR粒子滤波算法和IPF粒子滤波算法作为对照组.

初始化参数, 隐状态实际初始位置x₀ = 0.1, 初始化粒子集x_s~N(0.3, 8), 粒子集大小N = 100, 估计步数K= 75, IPF算法预设交叉权值α_c = 0.1, 突变权值α_m= 0.5.

实验分为非线性动态模型隐状态跟踪和非线性状态突变动态模型隐状态跟踪两部分.

4.1 非线性动态模型隐状态跟踪

对3种算法跟踪隐状态的性能进行仿真, 得到如下结果.

本文算法解决的核心问题是粒子贫化.针对该问题, 提出衡量参数N_S, 该参数表示在一步估计中, 粒子集中的粒子种类数.因为本文实验预设粒子集的大小为100, 所以取N_S ∈ [0, 100].对每步跟踪粒子多样性的仿真结果如图 2所示.

由图 2可知:对比SIR算法, 本文算法的粒子多样性也有明显提升; 对比IPF算法, 本文算法的粒子也有提升, 并且使粒子种数以更小的方差稳定在较高的数值, 波动明显小于IPF算法.实验结果验证了本文算法具有提升和自适应保持粒子多样性的作用.

图 3为3种算法对目标隐状态的跟踪轨迹.可以看出, 本文算法对目标状态轨迹的跟踪精度最高.

图 4为3种算法各步的RMSE.对3种算法进行500次蒙特卡洛仿真, 得到所有采样点的RMSE均值如表 1所示.

综合图 4和表 1可知, 本文算法的跟踪精度最高.

4.2 非线性动态突变模型隐状态跟踪

基于式(20)和(21), 构建含状态突变的非线性动态模型, 如下所示:

状态转移方程为

(22)

状态观测方程为

(23)

其中:转移噪声为N_t~N(0, 1), 观测噪声为R_t~N(0, 2).减小噪声可以提高跟踪轨迹的重合度, 从而更清晰地反映状态突变处的变化.令k=15时出现突变, 得到仿真结果如图 5所示.

由图 5可知, 在波动出现后, 本文算法的粒子种数在最短步数内恢复到较高的数值上, 这验证了本文算法的自适应性和更高的鲁棒性.

图 6为3种算法对突变模型的隐状态跟踪轨迹图.可以看出: SIR算法在状态突变时误差较大, 多步后才恢复跟踪; IPF算法在突变处也出现较大跟踪误差, 调整速度较SIR算法有提升; 本文算法在动态模型突变时误差最小, 调整速度最快.

图 7为3种算法各步的RMSE, 可见在状态突变处, 本文算法的RMSE明显小于对照算法.对3种算法进行500次蒙特卡洛仿真, 得到所有采样点的RMSE均值如表 2所示.

综合图 7和表 2可知, 本文算法对状态突变模型的跟踪精度最高, 稳定性最好.

综上可知, 当状态出现突变时, 3种算法的粒子多样性均会下降.粒子多样性与跟踪精度同时下降, 表明系统状态的突变会加重粒子的贫化, 且粒子多样性直接影响算法的精度.实验验证了本文算法可自适应地保持粒子多样性, 对状态突变模型隐状态跟踪时, 具有更高的精度和鲁棒性.

5 结论

本文针对粒子滤波算法在跟踪非线性状态突变系统隐状态时, 因粒子贫化出现的估计精度下降问题, 提出了一种基于Student's t分布的自适应智能粒子滤波算法.该算法通过自适应地生成大小权值粒子集, 执行自适应加权交叉和突变操作, 提升粒子多样性, 从而提升粒子滤波算法的精度和鲁棒性.仿真结果表明:本文算法在一般非线性模型下, 可有效提升粒子多样性, 提升算法精度; 在状态突变模型下, 可实现对突变状态的自适应调整, 使粒子多样性保持稳定, 从而获得更高的精度和鲁棒性.