基于粒子滤波的被动传感器多目标跟踪方法

摘要

本发明公开了一种基于粒子滤波的被动传感器多目标跟踪方法，它属于制导技术领域，主要解决现有多目标跟踪方法跟踪易发散，目标状态估计不准确的问题。该方法通过粒子群优化及样本混合采样算法对多目标的样本分布进行优化，并结合联合概率数据关联算法对多目标进行跟踪。首先，利用粒子群优化算法优化多目标联合样本的分布，使其聚集到真实目标出现概率较大的高似然区域；其次，利用这些样本计算目标与观测之间的关联概率及目标后验概率分布；最后，在重采样过程中，将联合样本权值按各目标样本似然分解到对应的目标样本中去，各目标按分解后的权值独立重采样，进一步优化目标样本的分布，提高目标跟踪的精度。

说明书

技术领域

本发明属于制导技术领域，涉及目标跟踪。具体地说是一种基于粒子群优化和序贯蒙特卡罗的被动传感器多目标跟踪方法，可用于红外制导等系统。

背景技术

多目标跟踪中，由于目标漏检和杂波的影响，传感器所得量测和目标之间的关联存在不确定性，并且在被动条件下所测的角度信息是目标状态的非线性函数，因此要准确估计目标状态实现目标跟踪，就需要解决量测和目标的数据关联以及非线性滤波两个问题。

传统的多目标跟踪方法包括最近邻法NN，联合概率数据关联JPDA，多假设跟踪MHT算法，其中最近邻法是将离目标状态最近的量测直接与目标相关联，当测量精度较高时，跟踪性能较好，当测量精度下降时，其跟踪性能也将严重下降；多假设跟踪则是穷举目标与量测之间所有可能的关联事件，并逐步按时间扩展，其缺点是计算时间将随目标数和量测数成指数增长；JPDA是到现在为止解决数据关联最有效的方法之一，它给每对目标和量测的关联赋予一定的概率，然后结合贝叶斯准则通过预测和更新两个步骤完成目标后验概率和状态的估计。

基于序贯蒙特卡罗的算法SMC是近年发展起来的非线性滤波方法，有学者将JPDA与SMC相结合，用以解决多目标跟踪问题，利用一定数量的样本和对应的权值来拟合运动目标的后验概率分布，从理论上讲，当所取样本数趋于无穷多时，SMC可以拟合任意概率分布。但是由于实际应用中，考虑到跟踪精度和实时性的综合要求，样本数通常有限，其在采样和重采样的过程中会出现样本贫化的现象，使得样本失去多样性，状态估计不稳定，导致跟踪发散。

发明内容

针对上述问题，本发明提出一种基于粒子滤波的被动传感器多目标跟踪方法，以保持样本的多样性，提高目标的跟踪精度。

实现本发明的技术关键是：利用粒子群优化算法优化多目标联合样本的分布，使其聚集到各目标状态的高似然区域，即真实目标出现概率较大的区域，这样用以滤波的样本将具有丰富的多样性且每个样本重要性得以提高；利用这些联合样本计算目标与量测之间的关联概率及目标滤波分布，并且在重采样过程中，不再按照多目标串联产生的联合样本的权值采样，而是将联合样本权值按各目标样本似然分解到对应的目标样本中去，进一步优化目标样本的分布，提高目标跟踪的精度，具体实现步骤包括如下：

(1)根据各目标的初始分布抽取目标样本，构造联合样本：

${\left\{x_0^n\right\}}_{n=1}^N={\{x_0^{n,1},L,x_0^{n,i},L,x_0^{n,c}\}}_{n=1}^N;$

其中，n表示联合样本序号，i表示目标序号，N表示联合样本个数，c表示目标个数，表示0时刻第n个联合样本中目标i的样本，每个联合样本的初始权值取为

(2)计算t时刻的预测联合样本：

i∈[1，c]，n∈[1，N]，t≥1，其中，是t时刻第n个联合样本中目标i的样本；

(3)按如下步骤优化粒子群：

(3a)将t时刻的预测联合样本中的各目标样本作为粒子群优化的初始样本为目标样本赋予初始速度：

(3b)计算t时刻目标样本对传感器量测的似然，表示为其中，k＝1，L，m是粒子群优化迭代序号，m≥5为设定的总的粒子群优化迭代次数；

(3c)根据第1到第k次迭代中各个目标样本的似然，找出目标i中各个样本的个体最优解

(3d)根据第i个目标中所有样本的似然，找出该目标所有样本中的全局最优解

(3e)利用粒子群优化算法中的更新方程求得目标样本在第k+1次迭代中的位置和速度

(3f)重复步骤(3b)～(3e)m次，得到粒子群优化后的联合样本：

${\{x_t^{n,1},L,x_t^{n,i},L,x_t^{n,c}\}}_{n=1}^N={\{x_{n,t}^{m,1},L{x}_{n,t}^{m,i}L,x_{n,t}^{m,c}\}}_{n=1}^N,$

其中，为优化后的目标样本；

(4)按如下步骤对联合样本权值更新及归一化：

(4a)根据优化后的目标样本所对应的量测值，计算目标i在t时刻量测的均值和方差选出满足的所有有效量测j∈[1，M_t]，其中，y_t为被动传感器获得的量测，ε＝9.21为设定的门限值，M_t为t时刻所有有效量测的个数；

(4b)列举出有效量测与目标i的关联事件φ_i，j；

(4c)计算有效量测与目标i基于样本形式的关联似然由目标运动的马尔科夫性和贝叶斯准则，计算第n个联合样本中边缘关联事件φ_i，j的概率：p(φ_i，j|Y^t)ⁿ，其中，Y^t表示从第1到第t时刻所有有效量测的集合；

(4d)求第n个联合样本的所有关联事件的概率的和，得到第n个联合样本的权值并对其归一化，得到归一化权值

(5)由联合样本及其对应的权值通过对联合样本加权求和估计出各目标状态，作为结果输出，并同时执行步骤(6)；

(6)按如下步骤对联合样本权值分解及重采样：

(6a)将第n个联合样本的归一化权值写成c个目标样本权值求和的形式：

${\bar{w}}_t^n={\bar{w}}_t^{n,1}+L{\bar{w}}_t^{n,i}L+{\bar{w}}_t^{n,c},$

其中第i个目标样本的权值通过第i个目标样本的似然计算获得；

(6b)从N个联合样本权值中，各取出第i个目标样本的权值构成根据这些权值，采样N个新的样本其中，样本其对应的权值为分别是t时刻目标i重采样前的第l个样本及其对应的权值；

(7)重复步骤(2)，继续跟踪目标。

本发明具有以下优点：

(1)本发明采样粒子群优化算法改善了目标样本的分布情况，使得目标样本向目标出现概率大的高似然区域聚集，且每个样本的重要性得以提高，在目标样本较少的条件下即可达到较高的跟踪精度；

(2)本发明考虑了目标相近存在相互影响及耦合的情况，对多目标样本进行混合采样，即在目标样本的重采样阶段，先将联合样本权值按各目标似然分解到对应的目标样本中去，再对每个目标按分解后的权值独立重采样，使各目标中大权值样本得以复制，小权值样本得以抑制，进一步优化了目标样本分布，提高跟踪精度。

附图说明

图1是本发明的整体流程图；

图2是本发明使用的粒子群优化粒子速度和位置更新示意图；

图3是用本发明进行一次目标跟踪的效果图；

图4是用本发明进行位置跟踪的均方根误差图。

具体实施方式

一、基础理论介绍

1.系统方程

笛卡尔坐标系下，系统状态取x，y方向的位置和速度，可建立如下的非线性动态系统模型：

$x_{t+1}^i={\mathrm{Fx}}_t^i+{\mathrm{Gv}}_t^i---1)$

$y_t=h\left(x_t^i\right){+e}_t---2)$

其中，i＝1，L，c表示目标的序号，c表示总的目标数目，分别表示目标i在x方向和y方向上的坐标，分别表示目标i在x方向和y方向上的速度，下标t∈N表示时间，状态噪声服从方差为的零均值高斯分布，F，G分别为状态转移矩阵和输入矩阵，h为非线性函数，量测噪声e_t服从方差为R的零均值高斯分布，与e_t相互独立，y_t为传感器的量测值。

本发明中假设被动传感器只能观测目标的方位角信息，因此h定义如下：

$h\left(x_t^i\right)@a\tan \frac{y_t^i-y_o}{x_t^i-x_o}---3)$

其中，x_o，y_o为传感器的位置。

2.粒子群优化

设在一个D维的搜索空间中有N个粒子组成的种群X＝{x₁，L x_n L，x_N}，其中，第n∈[1，N]个粒子的位置和速度分别为x_n＝(x_n1，x_n2，L x_nD)和v_n＝(v_n1，v_n2，L，v_nD)，且其位置的最优解为s_n＝(s_n1，s_n2，L，s_nD)，而整个种群位置的最优解为g＝(g₁，g₂，L g_D)，则第k次粒子群优化迭代中第n个粒子位置和速度的更新如下：

$v_{\mathrm{nd}}^{k+1}=v_{\mathrm{nd}}^k+c_1\zeta (s_{\mathrm{nd}}^k-x_{\mathrm{nd}}^k)+c_2\eta (g_d^k-x_{\mathrm{nd}}^k)---4)$

$x_{\mathrm{nd}}^{k+1}=x_{\mathrm{nd}}^k+v_{\mathrm{nd}}^{k+1}---5)$

其中，k＝1，L，m是粒子群优化迭代的序号，m是预先设定的粒子群优化迭代的总次数，d＝1，L，D表示粒子维数的序号，表示第n个粒子的位置中的第d维数据，表示第n个粒子的速度中的第d维数据，表示第n个粒子的位置中第d维数据的最优解，表示整个种群中所有粒子的位置中第d维数据的最优解，c₁和c₂是学习因子，其经典取值为(0，2)之间的正常数，ζ和η为(0，1)之间均匀分布的伪随机数；由于表示的是第n个粒子当前位置与其本身最优位置的差向量，所以c₁表征了第n个粒子向其自身最优位置搜索的能力；而则表示第n个粒子当前位置与整个种群中粒子最优位置的差向量，所以c₂表征了该粒子向整个种群最优位置搜索的能力；式4)中包含表示第n个粒子速度的更新还要取决于其迭代前的速度。基本粒子群优化算法粒子位置和速度的更新如附图2所示。

二、本发明基于粒子滤波的被动传感器多目标跟踪方法

参照图1，本发明的具体实施步骤包括如下：

步骤1.初始化目标样本

令初始时刻t＝0，根据目标i的初始分布抽取目标样本并串联构造联合样本i∈[1，c]，n∈[1，N]，N为抽取的样本数，c为目标数，其中，和分别表示第n个联合样本中第i个目标样本在x方向和y方向上的坐标，和分别表示第n个联合样本中第i个目标样本在x方向和y方向上的速度，第n个联合样本的初始权值取为

步骤2.计算t时刻的预测联合样本

根据t-1时刻的目标样本及状态方程1)计算t时刻的预测样本用这些预测样本构造t时刻的联合样本：t≥1，其中，表示t时刻第n个联合样本中的第i个目标样本。

步骤3.对预测的联合样本进行粒子群优化

(3.1)将t时刻的预测联合样本中的各目标样本作为粒子群优化的初始样本样本的初始速度为：

(3.2)计算样本对传感器量测y_t的似然表示为

$f_{n,t}^{k,i}=p\left(y_t\right|x_{n,t}^{k,i})=\frac{1}{\sqrt{2\pi \det \left(R\right)}}\exp (-\frac{1}{2}{(y_t-y_{n,t}^{k,i})}^{\prime }{\left(R\right)}^{-1}(y_t-y_{n,t}^{k,i}))---6)$

其中，k＝1，L，m为粒子群优化迭代序号，m≥5为设定的总的粒子群优化迭代次数，y_t为传感器获得的量测值，R为量测协方差阵，为由目标样本根据量测更新方程2)计算得到的量测值；

(3.3)找出在第1到第k次迭代中，第n个联合样本中第i个目标样本的的最小值，用表示，取其对应的样本作为第n个联合样本中第i个目标样本的个体最优解

(3.4)找出第i个目标中所有样本的的最小值，用表示，取其对应的样本作为第i个目标中所有样本的全局最优解

(3.5)结合和样本在第k次迭代中的速度根据粒子群优化的速度更新方程4)，更新样本在第k+1次迭代中的速度

$v_{n,t}^{k+1,i}=v_{n,t}^{k,i}+c_1\zeta (s_{n,t}^{k,i}-x_{n,t}^{k,i})+c_2\eta (g_{n,t}^{k,i}-x_{n,t}^{k,i})---7)$

(3.6)结合和根据粒子群优化的位置更新方程5)更新样本

$x_{n,t}^{k+1,i}=x_{n,t}^{k,i}+v_{n,t}^{k+1,i}---8)$

(3.7)重复步骤(3.2)～(3.6)m次，得到优化后的联合样本${\left\{x_t^n\right\}}_{n=1}^N={\{x_t^{n,1},L{x}_t^{n,i}L,x_t^{n,c}\}}_{n=1}^N={\{x_{n,t}^{m,1},L{x}_{n,t}^{m,i}L,x_{n,t}^{m,c}\}}_{n=1}^N.$

步骤4.联合样本权值更新并归一化

(4.1)根据样本计算目标i在t时刻量测的均值和协方差

${\bar{y}}_t^i=\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}{y}_t^{n,i}---9)$

${\sigma }_t^i=\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}(y_t^{n,i}-{\bar{y}}_t^i){(y_t^{n,i}-{\bar{y}}_t^i)}^{\prime }---10)$

其中，为目标样本对应的量测值；

(4.2)利用均值和协方差选出满足式11)条件的有效量测的一个集合

$y_t^j=\{y_t:{(y_t-{\bar{y}}_t^i)}^{\prime }{\left({\sigma }_t^i\right)}^{-1}(y_t-{\bar{y}}_t^i)\le ϵ\}---11)$

其中，j＝1，，L M_t，M_t表示有效量测的总数目，ε＝9.21为设定的门限值；

(4.3)列举出量测与目标i的关联事件φ_i，j；

(4.4)计算有效量测与目标i基于样本形式的关联似然

$p\left(y_t^j\right|x_t^{n,i})=\frac{1}{\sqrt{2\pi \det \left({\sigma }_t^i\right)}}\exp (-\frac{1}{2}{(y_t^j-y_t^{n,i})}^{\prime }{\left({\sigma }_t^i\right)}^{-1}(y_t^j-y_t^{n,i}))---12)$

(4.5)由目标运动的马尔科夫性和贝叶斯准则计算第n个联合样本中边缘关联事件φ_i，j的概率p(φ_i，j|Y^t)ⁿ，

$p{\left({\phi }_{i,j}\right|Y^t)}^n=\frac{1}{c}{P}_d^{c-c^0}{(1-P_d)}^{c^0}{P}_f^{M_t-(c-c^0)}\underset{(j,i)\in \phi }{\Pi }p\left(y_t^j\right|x_t^{n,i})---13)$

其中，P_f和P_d分别表示虚警概率和目标检测概率，c⁰是关联事件φ_i，j中未检测到的目标的个数，Y^t表示从第1到第t时刻所有有效量测的集合；

(4.6)求第n个联合样本的所有关联事件的概率的和，得到第n个联合样本的权值并对其归一化，得到归一化权值

$w_t^n=\mathrm{\Sigma p}{\left({\phi }_{i,j}\right|Y^t)}^n---14)$

${\bar{w}}_t^n=w_t^n/\underset{n=1}{\overset{N}{\Sigma }}{w}_t^n.---15)$

步骤5.目标状态估计

利用步骤3得到的联合样本及步骤4得到的联合样本权值按式16)估计目标状态，作为结果输出，并同时执行步骤6，

$x_t^i=\underset{n=1}{\overset{N}{\Sigma }}{x}_t^n·{\bar{w}}_t^n.---16)$

步骤6.对联合样本权值进行分解并重采样

(6.1)计算t时刻第n个联合样本中第i个目标样本的似然：

$l_t^{n,i}=\frac{1}{\sqrt{2\pi \det \left({\sigma }_t^i\right)}}\exp (-\frac{1}{2}{(y_t-y_t^{n,i})}^{\prime }{\left({\sigma }_t^i\right)}^{-1}(y_t-y_t^{n,i}))---17)$

(6.2)根据第n个联合样本中第i个目标样本的似然，计算第n个联合样本中第i个目标样本的权值：

${\bar{w}}_t^{n,i}=(l_t^{n,i}/\underset{n=1}{\overset{c}{\Sigma }}{l}_t^{n,i}){\bar{w}}_t^n---18)$

则第n个联合样本的归一化权值可以写成c个目标样本权值求和的形式：

${\bar{w}}_t^n={\bar{w}}_t^{n,1}+L{\bar{w}}_t^{n,i}L+{\bar{w}}_t^{n,c};---19)$

(6.3)从N个联合样本权值中，各取出第i个目标样本的权值构成根据这些权值，采样N个新的样本其中，样本其对应的权值为分别是t时刻目标i重采样前的第l个样本及其对应的权值。

步骤7.重复步骤2，继续跟踪目标。

本发明的效果可通过以下实验仿真进一步说明：

1.仿真条件及参数

仿真场景如图3所示，出现在仿真场景中的每一个目标的真实状态为x＝[x，v_x，y，v_y]′，x，y分别为每一个目标在笛卡尔坐标系x方向和y方向上的坐标，v_x，v_y分别为每一个目标在x方向和y方向上的速度。目标的状态方程和量测方程分别如式1)和2)所示，且每个目标均服从常速模型：

$F=\left(\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right),$$G=\left(\begin{array}{cc}{T}^2/2& 0\\ T& 0\\ 0& T^2/2\\ 0& T\end{array}\right)$

其中，T是采样时间间隔，传感器提供方位角信息，仿真参数如表1所示，

表1实验仿真参数

2.仿真内容及结果分析

在三个传感器条件下三个目标的纯方位角跟踪仿真实验，仿真实验中对比了本发明的跟踪方法与现有iMC-JPDA和jMC-JPDA两种多目标跟踪方法的位置的均方根误差RMSE及失跟率，仿真结果分别如图4和表2所示，其中：

图4(a)为样本数N＝30条件下，本发明方法与iMC-JPDA和jMC-JPDA的位置均方根误差的比较图；

图4(b)为样本数N＝50条件下，本发明方法与iMC-JPDA和jMC-JPDA的位置均方根误差的比较图；

图4(c)为样本数N＝80条件下，本发明方法与iMC-JPDA和jMC-JPDA的位置均方根误差的比较图；

图4(d)为样本数N＝100条件下，本发明方法与iMC-JPDA和jMC-JPDA的位置均方根误差的比较图；

从图4(a)～4(d)可以看出，随着样本数的增加，三种跟踪方法的RMSE均减小，但本发明的RMSE始终低于iMC-JPDA和jMC-JPDA方法的RMSE。

表2为本发明方法与现有iMC-JPDA和jMC-JPDA的失跟率比较，

表2本发明方法与iMC-JPDA和jMC-JPDA的失跟率比较

由表2可以看出，在相同样本数的条件下，本发明方法失跟率明显低于iMC-JPDA和jMC-JPDA跟踪方法，且当样本数N超过30时，本发明方法不再出现失跟的情况。