0 概述

无线传感器网络(Wireless Sensor Network, WSN)因其低成本、部署迅速、具有容错性和自组织性等优势, 被广泛应用于各种科学和商业用途的跟踪定位^[1-3], 如分布式目标跟踪、导航和控制以及智能交通系统中的交通监控等。然而, 由于WSN的复杂性以及目标与传感器节点的相对运动, 它们之间的信号强度和相对相位会产生随机变化, 从而使得系统的测量过程表现出明显的非高斯特性。具有非高斯特性的噪声被看作是一类闪烁噪声^[4], 其存在对实时目标跟踪有巨大影响。高斯线性模型通常使用卡尔曼滤波器获得最优估计^[5], 但在处理此类非高斯噪声背景下的目标跟踪问题时, 传统卡尔曼滤波器可能无法得到有效的状态估计。因此, 研究一种全新可靠的目标跟踪算法对解决WSN目标跟踪问题具有重要意义。

粒子滤波(Particle Filtering, PF)基于序列蒙特卡洛滤波方法, 其通过加权粒子集合近似状态的后验密度函数(Posterior Density Function, PDF)获得系统状态的最小方差估计, 是目前解决系统非线性非高斯问题的理想方案^[6]。由于标准PF的估计性能较差, 因此许多学者对其提出改进。文献[7]提出一种辅助粒子滤波(Auxiliary Particle Filtering, APF)方法, 将当前时刻测量纳入采样粒子的生成过程, 但在噪声极大的情况下APF的性能会急剧下降。文献[8]提出一种鲁棒正则化粒子滤波(Regularization Particle Filtering, RPF)方法, 通过混合本地RPF来减少多模态的损失。此外, 文献[9]将PF与蝙蝠算法结合来提升算法的整体性能, 文献[10]将一种自适应调整似然分布的PF应用到非线性非高斯条件下的导航系统中, 文献[11]则提出一种改进的无迹粒子滤波(Spherical Simplex Unscented Particle Divergence, SSUPF)方法, 并将其应用到深空环境的高精度天文导航作业中。

值得注意的是, 采样粒子的数量对PF性能有较大影响。文献[12]证明了过少的采样粒子会导致发散, 而过多的粒子会增加计算负担并限制导航的实时性。为提高PF的效率, 文献[13-14]提出一种基于Kullback-Leibler距离(Kullback-Leibler Distance, KLD)的自适应PF方法, 并将其用于移动机器人定位。文献[15]综述了许多自适应采样粒子数目的PF技术。上述研究表明, 通过控制粒子数量来保证估计质量同时降低计算成本是可行的。

本文提出一种基于KLD的自适应粒子滤波算法，用于非高斯噪声背景下的WSN目标跟踪。通过对重要性概率密度函数(Importance Probability Density Function, IPDF)进行迭代更新，缩小其与真实后验分布的差异，以提高滤波器跟踪精度。同时针对标准PF算法中粒子规模固定的问题，结合KLD采样自适应调整粒子数量，从而保证算法的计算效率。

1 标准PF算法 1.1 非线性滤波模型

滤波算法通过式(1)对运动目标建模。

$ {\mathit{\boldsymbol{x}}_{k + 1}} = f\left( {{\mathit{\boldsymbol{x}}_k},{\mathit{\boldsymbol{u}}_{k|k - 1}},{\mathit{\boldsymbol{w}}_k}} \right) $

(1)

其中, f(·)为已知的状态函数。根据上一时刻状态x_k、目标控制量u_k以及过程噪声w_k来更新状态x_k+1, u_k的误差并入w_k以限制状态维度。

系统测量模型如下:

$ {\mathit{\boldsymbol{z}}_k} = h\left( {{\mathit{\boldsymbol{x}}_k}} \right) + {\mathit{\boldsymbol{v}}_k} $

(2)

其中, z_k∈$\mathbb{R}$^z为k时刻的测量值, h(·)是测量函数, v_k是测量噪声, 包括地图误差、传感器误差和地形测量误差等。

为了应对WSN的复杂情况, 文献[16]提出一种有效的多传感器测量模型(Multi-Sensor Model, MSM)。

$ \mathit{\boldsymbol{z}}_k^l = {h_k}\left( {\mathit{\boldsymbol{x}}_k^s,\mathit{\boldsymbol{r}}_k^s} \right) + {\mathit{\boldsymbol{v}}_k} $

(3)

其中, x_k^s和r_k^s表示目标和传感器的真实位置, l=0, 1, …, L是传感器编号, v_k表示测量噪声, z_k^l是当前传感器的测量值。设Z_k=[z_k¹, z_k², …, z_k^l]^T表示传感器的联合测量, 实际测量值计算公式为:

$ \mathit{\boldsymbol{z}}_k^s = \sum\limits_{i = 1}^l {{\mathit{\boldsymbol{\mu }}_i}\mathit{\boldsymbol{Z}}_k^i} $

(4)

其中, μ_i是每个时刻已知的传感器权重值, 可参考文献[16]中的式(8)~式(10)。

1.2 非高斯噪声模型

现实生活中的高斯噪声分布十分少见, 而非高斯噪声却经常出现, 例如闪烁噪声。为解决该问题, 研究人员使用高斯混合来模拟状态空间系统中的非高斯噪声。然而, 高斯混合模型对极值极为敏感, 如果在噪声观测中存在极值可能会导致估计性能显著下降。文献[17]证明闪烁噪声可以由一个具有拖尾特性的噪声叠加到高斯噪声中来表示。因此, 本文将一个具有拖尾特性的噪声叠加到高斯噪声中, 用以模拟由强烈随机摆动而产生的非高斯噪声。由于在相同方差条件下, 拉普拉斯分布的拖尾特性强于高斯分布, 因此本文用拉普拉斯噪声来描述叠加噪声。

$ v\left( x \right) = \left( {1 - \tau } \right){v_g}\left( x \right) + \tau {v_l}\left( x \right) $

(5)

$ {v_g}(x) = \frac{1}{{\sqrt {2{\rm{ \mathit{ π} }}\sigma } }}{{\rm{e}}^{ - \frac{{{\mathit{\boldsymbol{x}}^2}}}{{2{\sigma ^2}}}}} $

(6)

$ {v_l}(x) = \frac{1}{{2\eta }}{{\rm{e}}^{ - \frac{{\left| x \right|}}{\eta }}} $

(7)

其中, τ＜＜1表示闪烁概率。

1.3 粒子滤波算法

PF是一种基于序列蒙特卡洛的滤波方法, 其使用一组采样粒子S_k:={x_kⁱ, ω_kⁱ}_i=1^N的加权样本集合来近似条件概率密度p(x_k|Z_k)。后验概率分布可表示为:

$ p\left( {{\mathit{\boldsymbol{x}}_k}|{\mathit{\boldsymbol{Z}}_k}} \right) \approx \sum\limits_{i = 1}^N {\omega _k^i} \delta \left( {{\mathit{\boldsymbol{x}}_k} - \mathit{\boldsymbol{x}}_k^i} \right) $

(8)

其中, δ(·)表示狄拉克函数, N表示粒子总数, ω_kⁱ表示第i个粒子的权重。

通常无法直接从后验概率密度p(x_k|Z_k)中采样。因此, 一般使用连续重要性采样方法从IPDF q(x_k|Z_k)间接采样。假设q(x_k|Z_k)可分解为:

$ q\left( {{\mathit{\boldsymbol{x}}_k}|{\mathit{\boldsymbol{Z}}_k}} \right) = q\left( {{\mathit{\boldsymbol{x}}_k}|{\mathit{\boldsymbol{x}}_{k - 1}},{\mathit{\boldsymbol{Z}}_k}} \right)q\left( {{\mathit{\boldsymbol{x}}_{k - 1}}|{\mathit{\boldsymbol{Z}}_{k - 1}}} \right) $

(9)

根据重要抽样原则^[14], 非正态化重要性权重更新方程如下:

$ \omega _k^i = \omega _{k - 1}^i\frac{{p\left( {{\mathit{\boldsymbol{z}}_k}|\mathit{\boldsymbol{x}}_k^i} \right)p\left( {\mathit{\boldsymbol{x}}_k^i|\mathit{\boldsymbol{x}}_{k - 1}^i} \right)}}{{q\left( {\mathit{\boldsymbol{x}}_k^i|\mathit{\boldsymbol{x}}_{k - 1}^i,{\mathit{\boldsymbol{Z}}_{1:k}}} \right)}} $

(10)

在实际应用中, 选取合适的IPDF q(x_k|Z_k)十分重要。文献[18]证明了最优IPDF是后验密度函数p(x_k|x_k－1, Z_k), 因此, q(x_k|Z_k)的选取应尽量接近后验概率分布。

序贯重要性采样(Sequential Importance Sampling, SIS)算法将预测概率密度函数p(x_k|Z_k)作为IPDF, 即:

$ q\left( {{\mathit{\boldsymbol{x}}_k}|{\mathit{\boldsymbol{x}}_{k - 1}},{\mathit{\boldsymbol{Z}}_k}} \right) = p\left( {{\mathit{\boldsymbol{x}}_k}|{\mathit{\boldsymbol{x}}_{k - 1}}} \right) $

(11)

将式(11)代入式(10), 权重的递归更新表达式为:

$ \omega _k^i \propto \omega _{k - 1}^ip\left( {{\mathit{\boldsymbol{z}}_k}|\mathit{\boldsymbol{x}}_k^i} \right) $

(12)

然而, SIS算法在经过多次迭代后, 大量粒子的权值会衰减。为了解决粒子退化问题, 文献[19-20]引入重采样步骤, 提出采样重要性重采样(Sampling Importance Resampling, SIR)算法, 将具有高权重的粒子重新采样, 并丢弃具有低权重的粒子, 再通过状态空间模型和测量模型生成新的粒子集合{$\mathit{\boldsymbol{\widetilde x}}_k^j,\widetilde \omega _k^j$ }_j=1^N_r。

2 基于KLD采样的自适应粒子滤波算法

目前, 多数PF在整个滤波过程中手动选择固定数量的采样粒子。然而在实际过程中, 目标状态的概率分布可能会随时间急剧变化, 从而导致常规PF效率较低。此外, 在滤波过程中, PF选择的IPDF通常只包含系统状态信息, 而忽略了当前的测量信息, 也会降低滤波结果的可靠性。

文献[10]提出一种基于KLD采样的自适应PF算法, 根据采样粒子近似分布与真实概率分布之间的KLD来动态调整粒子数量。该方法默认采样粒子从真实后验分布中选取, 忽略了重要性概率密度分布与真实分布之间的差异。

本文提出一种自适应PF算法, 在必要时通过迭代更新IPDF, 引入模拟退火算法和自适应退火参数, 使得算法中的粒子可以快速接近高似然区域, 从而大幅提高算法采样效率。在此基础上, 将上述自适应算法与KLD采样结合, 进而得到一种能够动态调整粒子数量的自适应PF算法.

2.1 KLD采样

KLD反映了2个概率分布p(x)和q(x)之间的差异, 可以表示为:

$ K(p,q) = \sum\limits_x p (x)\ln \frac{{p(x)}}{{q(x)}} $

(13)

假设真实概率分布是离散化的, 包含m个不同子空间b, 真实概率分布和基于采样的最大似然估计之间的KLD小于指定误差界限ε的概率为1－ρ, 即:

$ p(K(p,\hat p) < \varepsilon ) = 1 - \rho $

(14)

根据文献[9]可导出满足该条件的必要最小数量粒子为:

$ n = \frac{{m - 1}}{{2\varepsilon }}{\left[ {1 - \frac{2}{{9(m - 1)}} + \sqrt {\frac{2}{{9(m - 1)}}} {z_{1 - \rho }}} \right]^3} $

(15)

其中, z_1－ρ是标准正态分布1－ρ分位数的上界。

在常规KLD-PF中, 由于子空间尺寸固定, 若在整个滤波过程中均使用尺寸过小的子空间, 易在早期采样阶段引起粒子激增, 并产生较大的初始误差, 同时大幅提高算法运算量。为此, 文献[21]提出一种自适应调节子空间尺寸的方案, 即在状态空间较大时选取较大的子空间, 反之选取较小的子空间。因此, 算法能够在全局范围动态调整粒子数量。子空间尺寸m的计算公式为:

$ m = \lambda \sqrt {prod\left( {ceil\left( {\max \left( {\mathit{\boldsymbol{x}}_k^i} \right) - \min \left( {\mathit{\boldsymbol{x}}_k^i} \right)} \right)/\left( {2\varepsilon {n_{\max }}} \right)} \right.} $

(16)

其中, prod(·)表示向量各维元素的乘积, n_max是整个滤波过程允许的最大粒子数量, ε表示误差界限, λ是一个可调节的参数。

2.2 迭代自适应粒子滤波算法

为了在滤波过程中得到有效的IPDF并实时确定粒子规模, 本文提出一种基于迭代搜索得到全局范围最优粒子集合的算法, 通过对目标模型及其估计之间的输出观测误差进行判别, 以确定在当前时刻是否需要进行更新迭代。该算法的基本思想是在PF过程中通过迭代将采样粒子集合推向真实状态空间。

假设在k－1时刻的粒子集合S_k－1:={x_k－1ⁱ, ω_k－1ⁱ}_i=1^N已知, 首先可以从p(x_k|x_k－1)中得到k时刻状态x_k, 然后通过理想测量方程y_k=h(x_k)将其转换到测量空间。因此, 后验概率密度函数p(x_k|x_k－1, Z_k)可以由一个高斯分布来近似N(y_k, P_k), 其中, ${\mathit{\boldsymbol{\overline y}} _k} = \sum\limits_{i = 1}^N {\omega _{k - 1}^i} \mathit{\boldsymbol{y}}_k^i,{\mathit{\boldsymbol{P}}_k} = \sum\limits_{i = 1}^N {\omega _{k - 1}^i} \left( {\mathit{\boldsymbol{y}}_k^i - {{\mathit{\boldsymbol{\overline y}} }_k}} \right){\left( {\mathit{\boldsymbol{y}}_k^i - {{\mathit{\boldsymbol{\overline y}} }_k}} \right)^{\rm{T}}}$。预设定一个阈值ξ, 若(z_k－y_k)^TP_k^－1(z_k－y_k)＜ξ, 说明k－1时刻的采样粒子集有效, 否则根据式(17)更新IPDF(假设u=0)。

$ {\mathit{\boldsymbol{x}}_k} = f\left( {{\mathit{\boldsymbol{x}}_{k - 1}}} \right) + {\Delta _k} + {\mathit{\boldsymbol{w}}_{k - 1}} $

(17)

其中, Δ_{j, k}是一个修正补偿函数, 通过式(18)迭代计算。

$ {\Delta _{j,k}} = \left\{ {\begin{array}{*{20}{l}} {0,j = 0}\\ {{{\mathit{\boldsymbol{\hat x}}}_{k,j - 1}} - f\left( {{{\mathit{\boldsymbol{\hat x}}}_{k - 1}}} \right),j = 1,2, \cdots ,J} \end{array}} \right. $

(18)

其中, ${\mathit{\boldsymbol{\widehat x}}_{k,j - 1}}$是迭代j－1次的估计。将式(18)代入式(17)即可得到新的IPDFq_j(x_k|x_k－1, z_k)=p_w(x_k－f(x_k－1)－Δ_{j, k})。用新得到的IPDF代替标准PF算法中的p(x_k|x_k－1ⁱ)即可得到一种基于迭代过程的自适应粒子滤波算法AI-PF, 从而得到新的状态估计${\mathit{\boldsymbol{\widehat x}}_{k,j}}$和修正补偿函数Δ_{j+1, k}。基于上述过程进行迭代, 直到满足预定条件j=J。然而直接进行上述算法可能会导致无法更新状态方程, 这是因为直接从q_j(x_k|x_k－1, z_k)采样依然会引起粒子退化。有一种方案是向粒子权重添加一个小且正的参数δ(通常δ=10^－20), 但这会导致修正补偿函数Δ_{j, k}≈0。

本文在上述过程中引入模拟退火算法。假设R表示测量噪声的方差, 在迭代步数j很小时, 为其引入一个退火参数α_j, 使得方差变为α_jR。其中α_j≥1且随着迭代进行逐渐趋近于1。考虑将α_j的求解转化为以下优化问题:{x_{k, j}ⁱ}_i=1^N表示第j次迭代得到的粒子集合, 将其转换到测量空间得到{y_{k, j}ⁱ}_i=1^N, 并通过矩匹配得到测量噪声(如果是非高斯噪声)的一个高斯近似分布。因此, 粒子权重可表示为:

$ \omega _j^i \propto \exp \left[ {{{\left( {\mathit{\boldsymbol{y}}_{j,k}^i - {\mathit{\boldsymbol{z}}_k}} \right)}^{\rm{T}}}{{\left( {{\alpha _k}\mathit{\boldsymbol{R}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{y}}_{j,k}^i - {\mathit{\boldsymbol{z}}_k}} \right)} \right] = \exp \left( { - {d^i}{\beta _j}} \right) $

(19)

其中, dⁱ=(y_{j, k}ⁱ－z_k)^TR^－1(y_{j, k}ⁱ－z_k), β_j=1/α_j。

为了快速将粒子集合推向高似然区域, 应使第θ个粒子的权重尽可能大, 其中θ=$\mathop {{\rm{argmin}}}\limits_i $ (y_{j, k}ⁱ－z_k)^T·R^－1(y_{j, k}ⁱ－z_k), 等效于求式(20)的最大值。

$ J\left( {{\beta _j}} \right) = \frac{{\exp \left( { - {d^\theta }{\beta _j}} \right) + \delta }}{{\sum\limits_{i = 1}^N {\left[ {\exp \left( { - {d^i}{\beta _j}} \right) + \delta } \right]} }} $

(20)

由于J(β_j)的分母是$N\left[ {\exp ( - \overline d {\beta _j}) + \delta } \right]$, 其中$\overline d = {({\mathit{\boldsymbol{\overline y}} _{j,k}} - {\mathit{\boldsymbol{z}}_k})^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}({\mathit{\boldsymbol{\overline y}} _{j,k}} - {\mathit{\boldsymbol{z}}_k})$, 忽略常数N, 因此可以得到等价性能指标:

$ J\left( {{\beta _j}} \right) = \frac{{\exp \left( { - {d^\theta }{\beta _j}} \right) + \delta }}{{\exp \left( { - \bar d{\beta _j}} \right) + \delta }} $

(21)

式(21)对β_j求导可得:

$ \begin{array}{l} \bar d\exp \left( { - \bar d{\beta _j}} \right)\left[ {\exp \left( { - {d^\theta }{\beta _j}} \right) + \delta } \right] = \\ \;\;\;\;\;\;\;\;{d^\theta }\exp \left( { - {d^\theta }{\beta _j}} \right)\left[ {\exp \left( { - \bar d{\beta _j}} \right) + \delta } \right] \end{array} $

(22)

通过式(22)和不等式d>d^θ可以证明:

$ {{\rm{d}}^2}J\left( {{\beta _j}} \right)/{\rm{d}}\beta _j^2 < 0 $

考虑退火参数α_j≥1, 可由α_j=max{1, 1/β_j}计算得到α_j, 其中β_j是式(22)的解。

根据以上过程, 可以得到一个尽可能接近似然区域的IPDF。将上述方法和KLD采样结合, 在避免真实后验概率密度分布与IPDF之间差异过大的同时, 在KLD重采样过程中采用式(15)计算子空间数量。由文献[21]可知, 过大的子空间尺寸会导致过小的粒子集合, 从而导致较大的跟踪误差; 过小的子空间尺寸则会引发早期粒子激增。因此, 为了保证必要的跟踪性能, 在每次计算式(15)同时, 需为m设置上界和下界[m_min, m_max]。

综合上述改进, 本文提出一种修正的自适应迭代粒子滤波算法MAI-PF, 利用修正补偿函数提升算法性能的同时, 借助KLD采样减弱模拟退火算法对计算效率的负面影响。为了得到更加接近后验概率密度的IPDF, 算法在满足预设阈值条件时迭代更新系统方程, 即引入式(18)的修正补偿函数。算法的每一次迭代过程都是将粒子推向高似然区域的过程。由于似然函数的形状是尖锐的, 从q_j(x_k|x_k－1, z_k)中抽取的所有粒子都可能具有可忽略的权重, 导致迭代过程无法进行, 因此本文引入模拟退火算法和自适应退火参数α_j, 将当前测量结合到采样过程, 使得似然函数更加平滑, 从而大幅提高采样效率。自适应退火参数α_j可以通过计算式(22)得到。KLD可以衡量2个分布之间的差异, 满足PF对后验概率分布估计的应用情景。因此, 结合KLD采样可以降低粒子规模, 提高计算效率。通过式(15)计算子空间尺寸, 并为其设置必要的上下界, 进一步提升算法的自适应性。本文MAI-PF算法将KLD采样和模拟退火算法融入PF, 使其相互作用, 提高粒子多样性, 压缩粒子规模, 提高算法的精度与计算效率。MAI-PF算法的步骤为:

步骤1  输入预定阈值ξ、最小(最大)粒子数量N_min(N_max)、KLD误差界限ε、标准正态分布上分位数1－ρ以及子空间尺寸区间[m_min, m_max]。

步骤2  初始化粒子数n=0、子空间序号n_m=0、目标粒子集合S_k=Θ。

步骤3  在k－1时刻得到由初始IPDF生成的采样粒子集合S_k－1:={x_k－1ⁱ, ω_k－1ⁱ}_i=1^N, 同时根据预设阈值条件(z_k－y_k)^TP_k^－1(z_k－y_k)＜ξ判断是否需要引入补偿函数。若不需要, 则跳转步骤4, 否则跳转步骤3。

步骤4  更新IPDF。在迭代初期引入根据式(22)求得的退火参数α_j, 然后根据式(18)给出的补偿函数Δ_{j, k－1}来迭代更新IPDFq_j(x_k|x_k－1, z_k), 直到满足预设条件或迭代终止条件j=J, 并得到更新后的采样粒子集合S_{j, k－1}:={x_{j, k－1}ⁱ}_i=1^N。

步骤5  根据式(16)计算子空间尺寸并判断是否满足上下界条件m∈[m_min, m_max]。若m＜m_min, 令m=m_min; 若m>m_max, 令m=m_max。

步骤6  KLD重采样。从步骤3更新的粒子集合S_{j, k－1}中抽取粒子同时更新系统状态, 计算粒子权重, 并将新粒子添加到集合S_k=S_k∪{x_{k, j}ⁿ, ω_{k, j}ⁿ}。更新粒子数n=n+1, 判断新粒子是否落在空的子空间b_nm中。判断结果为是, 更新子空间数量n_m=n_m+1, 并将子空间b_nm标记为非空子空间。根据式(15)计算粒子数:

$ {n_{{\rm{KL}}}} = \frac{{{n_m} - 1}}{{2\varepsilon }}{\left[ {1 - \frac{2}{{9\left( {{n_m} - 1} \right)}} + \sqrt {\frac{2}{{9\left( {{n_m} - 1} \right)}}} {z_{1 - \rho }}} \right]^3} $

若n＜n_KL或者n＜n_min, 则跳转步骤6, 否则跳转步骤7。

步骤7  返回粒子集合S_k, 更新时刻k=k+1, 循环步骤2~步骤7。

3 仿真实验与结果分析 3.1 仿真条件与参数

本文实验建立跟踪目标与观测载机之间的运动模型如下:

$ {\mathit{\boldsymbol{x}}_{k + 1}} = {\mathit{\boldsymbol{F}}_k}{\mathit{\boldsymbol{x}}_k} + {G_k}{\mathit{\boldsymbol{w}}_k} $

(23)

其中, x=[r_x, v_x, r_y, v_y]^T表示运动目标的实时状态(r_x和r_y的单位为m, v_x和v_y的单位是m/s), F_k表示状态转移矩阵, G_k表示过程噪声转移矩阵, 噪声w_k为协方差为Q的0均值高斯白噪声。

$ {\mathit{\boldsymbol{F}}_k} = \left[ {\begin{array}{*{20}{c}} 1&T&0&0\\ 0&1&0&0\\ 0&0&1&T\\ 0&0&0&1 \end{array}} \right],{\mathit{\boldsymbol{G}}_k} = \left[ {\begin{array}{*{20}{c}} {{T^2}/2}&0\\ T&0\\ 0&{{T^2}/2}\\ 0&T \end{array}} \right] $

(24)

其中, T为采样周期。

本文测量模型采用第1.1节的MSM。在时刻k, 第l个观测载机测量表示为:

$ \mathit{\boldsymbol{z}}_k^l = {h_k}\left( {{\mathit{\boldsymbol{x}}_k},\mathit{\boldsymbol{x}}_k^{r,l}} \right) + \mathit{\boldsymbol{v}}_k^l $

(25)

$ {h_k}\left( {{\mathit{\boldsymbol{x}}_k},\mathit{\boldsymbol{x}}_k^{r,l}} \right) = \left[ {\begin{array}{*{20}{c}} {\sqrt {{{\left( {{r_{x,k}} - r_{x,k}^{r,l}} \right)}^2}} + {{\left( {{r_{y,k}} - r_{y,k}^{r,l}} \right)}^2}}\\ {\arctan \frac{{{r_{x,k}}}}{{{r_{y,k}}}}} \end{array}} \right] $

(26)

其中, x_k^{r, l}表示第l个观测载机在k时刻的状态, v_k=(v_r, v_α)为闪烁噪声, v_r、v_α相互统计独立, 其联合概率密度函数可表示为式(27)。

$ p\left( {{\mathit{\boldsymbol{v}}_k}} \right) = \varepsilon \prod\limits_{i = 1}^2 L \left( {{\mathit{\boldsymbol{v}}_i};0,{\eta _i}} \right) + (1 - \varepsilon )\prod\limits_{i = 1}^2 N \left( {{\mathit{\boldsymbol{v}}_i};0,\sigma _{{\mathit{\boldsymbol{v}}_i}}^2} \right) $

(27)

假设目标在1 000 m×1 000 m的区域内运动, 观测载机的数量为3, 其状态分别为x₁^r=[50, 0, 0, 0]^T, x₂^r=[0, 0, 150, 0]^T, x₃^r=[0, 10, 500, -10]^T。考虑一种实际场景, 在跟踪的起始时刻算法输入的并非目标真实状态。令协方差P₀=[49² 1 49² 1], 测量噪声由标准差为5的高斯噪声和η=1的拉普拉斯噪声叠加形成, 闪烁系数ε=0.1, 采样时间T=1 s, 蒙特卡洛仿真运行次数M=500, 预设阈值ξ=0.5, 退火迭代的终止条件为J=3, 在KLD-PF中的KLD误差界限ε和分位数1－ρ分别为0.15和0.99, 固定子空间尺寸为5 m(x、y方向一致), 本文算法的子空间尺寸界限m∈[5 m, 50m], λ的值为1, 最小和最大粒子数量分别为100和2 500。基于上述条件对标准PF、KLD-PF以及本文提出的MAI-PF 3种算法进行对比分析。

3.2 跟踪精度比较

本文使用目标位置的均方根误差(Root Mean Squared Error, RMSE)来评估算法的跟踪精度。

$ RMS{E_k} = \sqrt {\frac{1}{M}\left[ {{{\left( {\hat r_{x,k}^m - r_{x,k}^m} \right)}^2} + {{\left( {\hat r_{y,k}^m - r_{y,k}^m} \right)}^2}} \right]} $

(28)

其中, $\widehat r_{x,k}^m$和$\widehat r_{y,k}^m$表示在k时刻处的第m次蒙特卡洛仿真运行的估计位置, r_{x, k}^m和r_{y, k}^m是在k时刻处的真实位置, M是蒙特卡洛运行次数。对3种算法分别进行模拟实验, 其目标轨迹跟踪曲线如图 1所示。

为了验证初始粒子对3种算法的影响, 本文设计了2个实验进行对比分析。

实验1  选取500个初始粒子初始化算法, 3种算法的测试结果如图 2所示。

从图 2(a)可以看出, 3种算法的RMSE值非常接近, 表明3种算法能收敛。从图 2(b)可以看出, 随着时间的推移, KLD-PF和MAI-PF算法的粒子数量逐渐减少并趋于稳定, 标准PF的粒子数量保持在500个。此外, KLD-PF算法中存在一次迭代的粒子数量骤增, 意味着下一次迭代的时间增加, 而传感器采样周期是固定的, 这可能导致滤波过程中重要信息丢失。

实验2  选取100个初始粒子初始化算法, 3种算法的测试结果如图 3所示。从图 3可以看出, 由于KLD-PF与MAI-PF算法能够实时调整采样粒子数量, 因此即使在初始粒子数量较小时, 依然能够保持较高的跟踪精度, 且MAI-PF的精度优于KLD-PF, 而标准PF的性能明显下降。此外, 结合图 2(b)和图 3(b)可以看出, 与KLD-PF相比, MAI-PF在整个滤波过程中均能将粒子数量保持在较为理想的区间内。

3.3 适应性比较

实验3  为了研究不同WSN环境条件对算法性能的影响, 本文在实验1条件下调整闪烁噪声的参数, 令τ^*=1.5τ, η^*=3η, 其他条件不变, 对3种算法分别进行实验, 结果如图 4所示。从图 4可以看出, 本文算法的跟踪精度优于KLD-PF和标准PF算法。由于MAI-PF是通过服从状态后验概率分布的粒子集合来进行状态估计, 对于非高斯噪声的干扰不敏感, 因此其曲线并没有较大的波动。而标准PF的RMSE曲线却变得更加曲折, 且跟踪精度变低。因此, MAI-PF较标准PF具有更好的适应性。此外, 本文引入不同下界参数对比算法性能。方便起见, 本文实验假设目标初始状态已知并位于原点, 其他条件与实验1相同, 实验结果如图 5所示。其中, 图 5(a)中下界均为5 m, 图 5(b)中上界均为50 m。

从图 5(a)可知, 随着上界增大, 算法使用的粒子规模会随之减小, 从而导致跟踪精度显著下降。从图 5(b)可知, 随着下界降低, 单次迭代后的粒子规模会变得非常大, 从而明显增加计算量, 甚至会超过传感器的采样周期, 影响该时刻测量值。因此, 从理论上来说, 大的子空间尺寸会减少运算量, 小的子空间尺寸则会提高跟踪精度, 但是为了保持算法的整体跟踪性能, 须对其设置上下界, 具体数值应参考目标运动地图先验信息与所需跟踪精度。

3.4 计算量比较

在不同初始粒子数的条件下, 标准PF、无KLD采样过程的自适应迭代粒子滤波(AI-PF)以及与KLD采样结合的修正滤波(MAI-PF)3种算法在一次蒙特卡洛仿真中的单次迭代运算时间如表 1所示。从表 1可以看出, 由于AI-PF和MAI-PF均引入模拟退火算法来优化算法性能, 因此这2种算法的运算时间均比标准PF有一定程度的增加。对比固定粒子规模的AI-PF, MAI-PF可以动态调整粒子数, 能够明显减少运算时间, 这证明本文算法在保持稳定滤波性能的同时, 可以有效提高计算效率。

4 结束语

本文提出用于WSN目标跟踪的MAI-PF算法。通过设计迭代算法更新IPDF，并与KLD采样相结合动态调节采样粒子数量，从而提高跟踪性能。仿真结果表明，该算法能够根据环境自适应调节每次迭代的粒子数量，在保证跟踪精度的同时，避免标准PF算法计算量大的问题，可满足实际工程要求。后续将优化本文算法，通过考虑传感器的位置误差，进一步提升其跟踪性能。