一种可变数目机动目标跟踪方法

摘要

本发明涉及一种以多模粒子滤波框架为基础，在粒子状态预测与更新步骤中，依据粒子存在变量进行预测粒子状态集的采样，考虑当前观测值与机动目标状态粒子的关联程度，利用模糊拍卖算法与粒子群优化理论解决观测集与机动目标状态采样粒子集之间的关联问题，并给出机动目标出现与消失的判定准则，实现粒子权值更新；利用序贯重要性采样理论对混合采样粒子集进行重采样，得到包含模型信息和状态信息并能逼近每个机动目标状态后验概率分布的粒子集；考虑粒子存在变量的影响，按照目标模型概率进行粒子状态融合得到机动目标局部状态后验估计值和均方误差；最后，对各传感器局部航迹信息进行加权融合，得到各机动目标的全局状态估计值，实现对机动目标数目及状态变化的准确估计。

说明书

技术领域

本发明属于机动目标跟踪领域，涉及一种新的目标联合检测与跟踪方法，具体涉及一种可变数目机动目标跟踪方法。

背景技术

现代防空系统需要借助于多传感器提供的互补观测信息进行探测工作，复杂的观测环境对机动目标跟踪算法也提出了越来越高的要求。传统的先检测后跟踪算法大多先对传感器每帧原始数据进行门限检测，然后利用超过门限值的测量数据进行跟踪处理，这种方法虽然在一定程度上减少了计算量，但是却损失了许多有用信息。为了提高目标的检测与跟踪性能，利用未经门限处理的多帧原始观测数据进行能量积累，以实现融合检测与跟踪的算法受到了很大的关注。典型的检测与跟踪算法包括基于Hough变换的算法，动态规划算法，极大似然估计，以及粒子滤波算法等，其中，基于粒子滤波实现的递归算法得到了广泛的改进与研究应用。

多模粒子滤波算法利用粒子滤波实现对各粒子加权求和得到状态估计，然后依据多模型中每个模型的模型概率对各模型的估计值加权求和，得到多模型粒子滤波器的全局估计值，该方法通常是在所采用的滤波器完成对状态估计的同时，进行实时系统模型估计，适用于机动目标跟踪。关于经典多模粒子滤波算法的改进方法大致涉及粒子重要性采样、重采样、模型粒子数以及滤波框架四类，它们虽然在一定程度上解决了粒子多样性及算法计算量问题，但是大多没有考虑实际多传感器观测环境下机动目标数目变化时的跟踪情况，并且忽视了对状态粒子集与当前观测值的关联判定，降低了目标跟踪性能。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种可变数目机动目标跟踪方法，

技术方案

一种可变数目机动目标跟踪方法，其特征在于：在S个多传感器的公共探测区域内，确定时刻k存在最大数目为Γ_max个的机动目标数目以及每个机动目标状态，具体步骤如下：

步骤1：对于传感器S_s观测下的机动目标状态进行随机采样，记机动目标t的状态为其初始状态分布为采样粒子γ对应于模型状态的模型概率为采样控制变量为u(γ)，其中u(γ)～U(0,1)，当满足时，该时刻第γ个粒子所对应的采样模型状态为以概率在模型对应的状态分布中抽取粒子相应粒子权重取得到机动目标t的初始状态粒子集 

其中：传感器编号为s＝1,2,…S，S为传感器最大数目，分别表示初始时刻机动目标t的采样粒子γ所对应的运动状态、采样模型状态及相应权值，t＝0,1,…,Γ_k，Γ_k＝0,1,…,Γ_max，Γ_max为时刻k可能存在的机动目标最大数目，k＝0,1,2…，γ＝1,2,…,Π，Π为采样粒子总数，M表示模型状态总数； 

步骤2、在k-1时刻，k＝1,2,…，记表征机动目标t混合状态后验分布的粒子集为 依据粒子存在变量进行粒子状态预测，其中，表示目标不存在，表示目标存在，且Π_e为转换概率；分三种情况进行处理：

状态(a)：对于持续存在的机动目标采样粒子，此时利用k-1时刻粒子模型状态通过模型状态转移矩阵Π_m，由求得当前k时刻粒子模型状态依据模型转移概率${\pi }_{\mathrm{\alpha \beta }}=P\{m_k^t=\beta |m_{k-1}^t=\alpha \},\alpha ,\beta =\mathrm{1,2},...,M,$k-1时刻对应模型状态α的模型概率可得时刻k机动目标t关于模型β的预测模型概率以预测模型概率为基础，根据不同模型状态对应的目标状态转移密度函数进行粒子状态预测，抽取Π_c个持续存在的机动目标采样粒子γ_c，形成粒子集

状态(b)：对于新生的机动目标采样粒子，此时新生机动目标的建议分布选取为初始状态分布，根据不同模型状态对应的状态分布抽取Π_b个新生的机动目标采样粒子γ_b，形成粒子集

状态(c)：对于消失的机动目标采样粒子，此时时刻k的粒子状态不存在；

步骤3：机动目标粒子集${\left\{{\tilde{x}}_k^{t,{\gamma }_c},{\tilde{e}}_k^{t,{\gamma }_c},{\tilde{m}}_k^{t,{\gamma }_c},{\tilde{w}}_k^{t,{\gamma }_c}\right\}}_{{\gamma }_c=1}^{{\Pi }_c}$和${\left\{{\tilde{x}}_k^{t,{\gamma }_b},{\tilde{e}}_k^{t,{\gamma }_b},{\tilde{m}}_k^{t,{\gamma }_b},{\tilde{w}}_k^{t,{\gamma }_b}\right\}}_{{\gamma }_b=1}^{{\Pi }_b}$记为 其中：为时刻k第η个机动目标状态采样粒子，i＝b,c，1≤η≤Γ_k，Γ_k为时刻k的机动目标数目；

对当前的观测值集合与机动目标状态采样粒子集合依据模糊拍卖算法和粒子群优化理论进行关联判定，并完成机动目标是否新生与消失的判断，为时刻k传感器S_s输出的第ζ个机动目标观测值，1≤ζ≤Ω_k，Ω_k为观测值数目；

具体如下： 

在k时刻，记第ζ个观测值与第η个机动目标状态采样粒子相关联的模糊代价为c_ζ,η，定义模糊隶属度为${\mu }_{{\Theta }_{\eta },k}=\exp \left\{{(z_{s,k}^{\zeta }-{\bar{z}}_{s,k}^{\eta })}^T{\Sigma }_{\eta ,k}^{-1}(z_{s,k}^{\zeta }-{\bar{z}}_{s,k}^{\eta })\right\},$其中， ${\bar{z}}_{s,k}^{\eta }=\frac{1}{{\Pi }_i}\underset{{\gamma }_i=1}{\overset{{\Pi }_i}{\Sigma }}{z}_{s,k}^{\eta ,{\gamma }_i},{\Sigma }_{\eta ,k}^{-1}=\frac{1}{{\Pi }_i-1}\underset{{\gamma }_i=1}{\overset{{\Pi }_i}{\Sigma }}(z_{s,k}^{\eta ,{\gamma }_i}-{\bar{z}}_{s,k}^{\eta }){(z_{s,k}^{\eta {,\gamma }_i}-{\bar{z}}_{s,k}^{\eta })}^T,$为机动目标状态采样粒子对应的预测观测值，模糊代价c_ζ,η可表示为对应的模糊隶属度观测值集合与机动目标状态采样粒子集合二维分配问题数学模型中的目标函数与约束条件为：

目标函数$F\left(a\right)=\min \underset{\zeta =1}{\overset{{\Omega }_k}{\Sigma }}\underset{\eta =1}{\overset{{\Gamma }_k}{\Sigma }}{c}_{\zeta ,n}·a_{\zeta ,\eta }$

约束条件$\left(\begin{array}{cc}\underset{\eta =1}{\overset{{\Gamma }_k}{\Sigma }}{a}_{\zeta ,\eta }\le 1& \forall 1\le \zeta \le {\Omega }_k\\ \underset{\zeta =1}{\overset{{\Omega }_k}{\Sigma }}{a}_{\zeta ,\eta }\le 1& \forall 1\le \eta \le {\Gamma }_k\\ a_{\zeta ,\eta }=\mathrm{0,1}& 1\le \zeta \le {\Omega }_k,1\le \eta \le {\Gamma }_k\end{array}\right)$

其中：a_ζ,η为二进制决策变量，a_ζ,η＝1表示第ζ个观测值与第η个机动目标状态采样粒子相关联，a_ζ,η＝0表示其他，该二维分配问题的最优解可表示为每行、每列最多只有一个1的Ω_k×Γ_k阶关联决策矩阵

定义粒子γ_pso的适应度值为其中，j_ζ,η为依据拍卖算法理论定义的竞拍价格，j_ζ,η＝λ₁g_ζ,η-λ₂v_ζ+ε，g_ζ,η为第ζ个观测值与第η个机动目标状态粒子集的关联价值，g_ζ,η＝max{c_ζ,η}-c_ζ,η，v_ζ为观测值的价值，λ₁为利益系数，λ₂为代价系数，ε为足够小的正数，设定粒子为个体当前最优位置，p_g为群体中最佳粒子位置，粒子的位置和速度向量维数为Ω_k×Π_pso，Ω_k为观测集维数，Π_pso为粒子群的粒子数，第γ_osp个粒子的位置表示一组观测值与机动目标状态采样粒子的关联性判断结果；

首先依据关联决策矩阵进行机动目标是否新生与消失的粗判断：若关联决策矩阵中某一行的元素值均近似为零，则认为该观测值可能来自于新生目标，或者来自于杂波；若关联决策矩阵中某一列的元素值均近似为零，则认为该目标可能消失；

然后依据目标在观测场景中的存在概率进行精判断：若P_e≥2/3，则时刻k新出现的观测值来自于新生目标t，此时在整个粒子集中增加该新生目标粒子集；若P_e≤1/3，则认为时刻k目标t消失，此时在整个粒子集中删除该消失目标粒子集；所述 $P_e=\frac{1}{{\Pi }_i}\underset{{\gamma }_i=1}{\overset{{\Pi }_i}{\Sigma }}{\tilde{e}}_k^{\eta {,\gamma }_i};$

步骤4：依据观测分辨单元内机动目标存在与不存在的似然比表达式 $L\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma }),L\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma })=\left(\begin{array}{cc}\frac{p\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma }=1)}{p\left(z_{s,k}\right|e_k^{t,\gamma }=0)}& e_k^{t,\gamma }=1\\ 1& e_k^{t,\gamma }=0\end{array}\right)$进行粒子权值计算，

其中：为机动目标状态采样粒子存在的似然函数，为机动目标状态采样粒子不存在的似然函数；

分三种情况进行：

状态(a)：对于持续存在的机动目标采样粒子，依据预测粒子状态时判定出的目标模型状态及相应的目标状态转移密度函数，得持续存在的机动目标采样粒子γ_c的未归一化权值为${\tilde{w}}_k^{t,{\gamma }_c}=w_{k-1}^{t,\gamma }L\left(z_{s,k}\right|x_k^{t,{\gamma }_c},e_k^{t,{\gamma }_c}=1),$归一化权值为$w_k^{t,{\gamma }_c}=\frac{{\tilde{w}}_k^{t,{\gamma }_c}}{\underset{{\gamma }_i=1}{\overset{{\Pi }_c}{\Sigma }}{\tilde{w}}_k^{t,{\gamma }_i}};$

状态(b)：对于新生的机动目标采样粒子，依据预测粒子状态时判定出的目标模型状态及相应的目标状态分布，得新生的机动目标采样粒子γ_b的非归一化权重为${\tilde{w}}_k^{t,{\gamma }_b}=\frac{L\left(z_{s,k}\right|x_k^{t,{\gamma }_b},{\tilde{e}}_k^{t,{\gamma }_b}=1)p\left(x_k^{t,{\gamma }_b}\right|{\tilde{e}}_k^{t,{\gamma }_b}=1,e_{k-1}^{t,\gamma }=0)}{{\Pi }_{b}q\left(x_k^{t,{\gamma }_b}\right|{\tilde{e}}_k^{t,{\gamma }_b}=1,e_{k-1}^{t,\gamma }=0,z_{s,k})},$归一化权值为$w_k^{t,{\gamma }_b}=\frac{{\tilde{w}}_k^{t,{\gamma }_b}}{\underset{{\gamma }_i=1}{\overset{{\Pi }_b}{\Sigma }}{\tilde{w}}_k^{t,{\gamma }_i}},$其中， $p\left(x_k^{t,{\gamma }_b}\right|{\tilde{e}}_k^{t,{\gamma }_b}=1,e_{k-1}^{t,\gamma }=0)$为新生机动目标采样粒子先验分布，$q\left(x_k^{t,{\gamma }_b}\right|{\tilde{e}}_k^{t,{\gamma }_b}=1,e_{k-1}^{t,\gamma }=0z_{s,k})$为新生机动目标采样粒子建议分布；

状态(c)：对于消失的机动目标采样粒子，粒子权值不存在；

将持续存在的机动目标采样粒子与新生的机动目标采样粒子的粒子集进行混合，记机动目标新生概率为P_b，消失概率为P_d，k-1时刻机动目标存在概率为可得非归一化混合概率${\tilde{M}}_b=P_b[1-{\hat{P}}_{k-1}]\underset{{\gamma }_b=1}{\overset{{\Pi }_b}{\Sigma }}{\tilde{w}}_k^{t,{\gamma }_b},{\tilde{M}}_c=[1-P_d]{\hat{P}}_{k-1}\underset{{\gamma }_c=1}{\overset{{\Pi }_c}{\Pi }}{\tilde{w}}_k^{t,{\gamma }_c},$以及k时刻机动目标的存在概率${\hat{P}}_k=\frac{{\tilde{M}}_b+{\tilde{M}}_c}{{\tilde{M}}_b+{\tilde{M}}_c+P_d{\hat{P}}_{k-1}+[1-P_b][1-{\hat{P}}_{k-1}]},$更新持续存在的机动目标采样粒子与新生的机动目标采样粒子权重可得合并Π_b+Π_c个粒子，得k时刻机动目标混合粒子集 其中i＝b,c，t＝0,1,…,Γ_k；

步骤5：对步骤4获得的混合粒子集中的Π_b+Π_c个粒子进行重采样，得到逼近机动目标t状态后验概率分布的粒子集其中，k＝1,2,…，t＝0,1,…Γ_max， $w_k^{t,\gamma }=\frac{1}{\Pi };$

步骤6：在时刻k，按照每个模型的模型概率进行粒子状态融合，得到机动目标t的局部状态后验估计值均方误差$P_{k/k}^t=\underset{m_k^{t,\gamma }=1}{\overset{M}{\Sigma }}\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}{e}_k^{t,\gamma }{w}_k^{t,\gamma }({\hat{x}}_{k/k}^{m_k^{t,\gamma },t}-{\hat{x}}_{k/k}^t){({\hat{x}}_{k/k}^{m_k^{t,\gamma },t}-{\hat{x}}_{k/k}^t)}^T,$其中：对应于模型${\hat{x}}_{k/k}^{m_k^{t,\gamma },t}=\frac{\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}{e}_k^{t,\gamma }{x}_k^{t,\gamma }{w}_k^{t,\gamma }}{\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}{e}_k^{t,\gamma }},\gamma =\mathrm{1,2},...,\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }},$为预测模型概率，为对应模型的采样粒子数，${\mu }_k^{m_k^{t,\gamma }}=\frac{{\mu }_{k/k-1}^{m_k^{t,\gamma }}{L}_k^{m_k^{t,\gamma }}}{\underset{m_k^{t,\gamma }=1}{\overset{M}{\Sigma }}{\mu }_{k/k-1}^{m_k^{t,\gamma }}{L}_k^{m_k^{t,\gamma }}},$为模型似然，$L_k^{m_k^{t,\gamma }}=\frac{\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}p\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma }=1)}{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}};$

依据目标存在的后验概率表达式P_e，得到时刻k机动目标数目Γ_k的估计值${\hat{\Gamma }}_k=\underset{t=1}{\overset{{\Gamma }_k}{\Sigma }}\left(\mathrm{Round}(P_e-\frac{1}{6})\right),$其中：$P_e\stackrel{\Delta }{=}P(e_k^{t,\gamma }=1|z_{s,k})=\frac{1}{\Pi }\underset{\gamma =1}{\overset{\Pi }{\Sigma }}{e}_k^{t,\gamma },$Round(x)表示求得距离x最近的整数；

在观测时间序列k＝1,2,…K内，对于传感器S_s，机动目标t的子航迹可表示为状态集合时刻k估计出的机动目标数目为

步骤7：重复步骤2至步骤6，求得各传感器观测下机动目标t的S个局部状态估计值及其相应的估计均方误差阵参照分布式融合结构， 得全局状态估计值全局状态估计均方误差观测场景中的目标数目估计值为各传感器局部估计值中的最大值，${\hat{\Gamma }}_{g_k}=\max \{{\hat{\Gamma }}_{1_k},{\hat{\Gamma }}_{2_k},...,{\hat{\Gamma }}_{s_k}\}.$

有益效果

本发明提出的一种可变数目机动目标跟踪方法，以多模粒子滤波框架为基础，在粒子状态预测与更新步骤中，依据粒子存在变量进行预测粒子状态集的采样，考虑当前观测值与机动目标状态粒子的关联程度，利用模糊拍卖算法与粒子群优化理论解决观测集与机动目标状态采样粒子集之间的关联问题，并给出机动目标出现与消失的判定准则，实现粒子权值更新；利用序贯重要性采样理论对混合采样粒子集进行重采样，得到包含模型信息和状态信息并能逼近每个机动目标状态后验概率分布的粒子集；考虑粒子存在变量的影响，按照目标模型概率进行粒子状态融合得到机动目标局部状态后验估计值和均方误差；最后，对各传感器局部航迹信息进行加权融合，得到各机动目标的全局状态估计值，实现对机动目标数目及状态变化的准确估计。

与传统多模粒子滤波算法相比，本发明在采样过程中使用了一定的先验信息和前一时刻的后验信息，通过对多模型混合采样得到机动目标状态后验分布粒子集，保证了粒子多样性；该方法引入了观测值与机动目标状态采样粒子之间的关联判定，给出机动目标出现与消失的判定准则，实现了对机动目标运动模型及状态变化的准确估计；同时，考虑粒子存在变量影响，按照每个模型的模型概率进行粒子状态融合得到机动目标状态估计值，提高了机动目标数目与状态估计的准确性，优化了目标跟踪性能。

附图说明

图1是本发明结构原理图；

表1是目标出现的时间段及相应运动状态；

图2是不同时刻各目标的实际位置

(a)不同时刻各目标在二维空间上的实际位置；(b)不同时刻各目标在一维空间上的实际位置；

图3是不同方法对目标1、2的运动模型概率估计比较结果

(a)本发明方法对目标1的运动模型概率估计；(b)方法一对目标1的运动模型概率估计；(c)方法二对目标1的运动模型概率估计；(d)本发明方法对目标2的运动模型概率估计；(e)方法一对目标2的运动模型概率估计；(f)方法二对目标2的运动模型概率估计；

图4是不同方法对目标1、2的位置估计误差比较结果

(a)不同方法对目标1的X方向位置状态估计；(b)不同方法对目标1的Y方向位置状态估计；(c)不同方法对目标2的X方向位置状态估计；(d)不同方法对目标2的Y方向位置状态估计；

图5是不同方法对目标数目估计的比较结果。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

用于实施的硬件环境是：因特尔酷睿2双核2.93G计算机、2.0GB内存、512M显卡，运行的软件环境是：Matlab R2012b，Windows 7。我们用Matlab R2012b软件实现了本发明提出的方法。

本发明具体实施如下：

步骤1：对机动目标状态和滤波初值进行初始化，具体为：为简化问题，在二维监视区域内，假设多传感器系统由一部雷达和一台红外探测器组成，各传感器布置在不同地点，不考虑传感器漏检现象，假设其采样时间同步，且已完成坐标转换与时间对准，设定所有目标在传感器观测的重合区域内机动。假设观测时间为150s，在整个观测期 间，共有三个目标出现和消失，目标在场景中出现的时间段及相应机动状态如表1所示。目标1的初始位置为[10500,8000]m，初始速度为[50,180]m/s；目标2的初始位置为[12500,15300]m，初始速度为[154,115]m/s；目标3的初始位置为[12100,9300]m，初始速度为[155,88]m/s；三个目标机动转弯均采用了两种转弯速率0.06rad/s与0.07rad/s。为了计算方便，设公共观测区域为[10000,27000]m×[4000,21000]m的观测空间，对于雷达传感器S₁，目标的距离单元均匀分布在[13000,30000]m，多普勒单元均匀分布在[-125,200]m/s，考虑只有一个方位单元的情况，分辨单元数目设为340×50×1，距离分辨率单元为50m，多普勒分辨率单元为6.5m/s，扫描间隔T₁＝1s。对于红外传感器S₂，灰度图像在两个方向上的距离分辨单元为100m，分辨单元数目为170×170，扫描间隔为T₂＝1s。对于雷达传感器得到的观测图像，假设三个目标所在的初始距离-多普勒分辨单元依次为(5,43)、(137,48)与(49,48)。对于红外传感器得到的观测图像，假设三个目标初始位置依次为(5,40)、(26,114)与(23,54)，速度依次为(0.5,1.8)、(1.53,1.11)与(1.65,0.92)分辨单元每秒。

表1

目标 出现时间段 出现时对应机动状态 转弯时间段 目标1 [1,100]、[111,150] 匀速与转弯、匀速 [51,90] 目标2 [31,120]、[123,140] 匀速与转弯、匀速 [71,100] 目标3 [21,130] 匀速与转弯 [61,110]

为验证本发明方法的有效性，使用其他两种滤波方法对场景中的机动目标进行跟踪，方法一为经典多模粒子滤波算法，方法二为一种改进的多模型粒子滤波算法。在各多模型方法中，均使用CV(Constant velocity)模型和CT(Constant turn)模型对目标进行估计。对于目标采样粒子对应的各初始模型概率均取为1/3，各目标模型概率的转移矩阵均取为Π_m＝[0.9,0.05,0.05；0.05,0.9,0.05；0.05,0.05,0.9]，目标初始存在概率取为0.5，目标新生与消失概率均取为0.05，粒子数均取为2000。新生目标粒子状态采样为所有分辨单元内的均匀采样，其中，粒子在X轴位置均匀分布在[10000,27000]m之间，速 度均匀分布在[-60,200]m/s之间；Y轴位置均匀分布在[4000,21000]m之间，速度均匀分布在[-200,200]m/s之间；目标强度均匀分布在[10,30]之间。设定关联优化算法中利益系数为0.8，竞价系数为0.2，种群规模为100，惯性权重为0.8，学习因子为2，迭代次数为50。由Matlab2012b仿真得出观测场景中机动目标的运动轨迹，如图2所示，其中，从图2(b)可以很直观的看出各目标在不同时刻的存在状态及相互关联程度。

步骤2：对于传感器S_s观测下的机动目标状态进行随机采样，记机动目标t的状态为其初始状态分布为采样粒子γ对应于模型状态的模型概率为采样

制变量为u(γ)，其中u(γ)～U(0,1)，当满足时，则认为该时刻第γ个粒子所对应的采样模型状态为以概率在模型对应的状态分布中抽取粒子相应粒子权重取得到机动目标t的初始粒子集 其中，传感器编号为s＝1,2,…S，S为传感器最大数目，分别表示初始时刻机动目标t的采样粒子γ所对应的运动状态、采样模型状态及相应权值，t＝0,1,…,Γ_k，Γ_k＝0,1,…,Γ_max，Γ_max为时刻k可能存在的机动目标最大数目，k＝0,1,2…，γ＝1,2,…,Π，Π为采样粒子总数，M表示模型状态总数。步骤3：在k-1时刻，k＝1,2,…，记表征机动目标t混合状态后验分布的粒子集为 依据粒子存在变量进行粒子状态预测，其中，表示目标不存在，表示目标存在，且Π_e为转换概率。分三种情况进行处理，具体描述如下：

状态(a)：对于持续存在的机动目标采样粒子，此时利用k-1时刻粒子模型状态通过模型状态转移矩阵Π_m，由求得当前k时刻粒子模型状态依据模型转移概率${\pi }_{\mathrm{\alpha \beta }}=P\{m_k^t=\beta |m_{k-1}^t=\alpha \},\alpha ,\beta =\mathrm{1,2},...,M,$k-1时刻对应模型状态α的模型概率可得时刻k机动目标t关于模型β的预测模型概率以预测模型概率为基础，根据不同模型状态对应的目标状态转移密度函数进行粒子状态预测，抽取Π_c个持续存在的机动目标采样粒子γ_c，得粒子集

状态(b)：对于新生的机动目标采样粒子，此时新生机动目标的建议分布选取为初始状态分布，根据不同模型状态对应的状态分布抽取Π_b个新生的机动目标采样粒子γ_b，得粒子集

状态(c)：对于消失的机动目标采样粒子，此时粒子状态不存在。

步骤4：将步骤3得到的机动目标粒子集记为其中，为时刻k第η个机动目标状态采样粒子，i＝b,c，1≤η≤Γ_k，Γ_k为时刻k的机动目标数目，对当前的观测值集合与机动目标状态采样粒子集合依据模糊拍卖算法和粒子群优化理论进行关联判定，并完成机动目标是否新生与消失的判断，为时刻k传感器S_s输出的第ζ个机动目标观测值，1≤ζ≤Ω_k，Ω_k为观测值数目。具体描述如下：

在k时刻，记第ζ个观测值与第η个机动目标状态采样粒子相关联的模糊代价为c_ζ,η，c_ζ,η的值越小说明第ζ个观测值与第η个机动目标状态采样粒子关联的可能性越大。定义模糊隶属度为${\mu }_{{\Theta }_{\eta },k}=\exp \left\{{(z_{s,k}^{\zeta }-{\bar{z}}_{s,k}^{\eta })}^T{\Sigma }_{\eta ,k}^{-1}(z_{s,k}^{\zeta }-{\bar{z}}_{s,k}^{\eta })\right\},$其中，${\bar{z}}_{s,k}^{\eta }=\frac{1}{{\Pi }_i}\underset{{\gamma }_i=1}{\overset{{\Pi }_i}{\Sigma }}{z}_{s,k}^{\eta ,{\gamma }_i},$${\bar{z}}_{s,k}^{\eta }=\frac{1}{{\Pi }_i}\underset{{\gamma }_i=1}{\overset{{\Pi }_i}{\Sigma }}{z}_{s,k}^{\eta ,{\gamma }_i},{\Sigma }_{\eta ,k}^{-1}=\frac{1}{{\Pi }_i-1}\underset{{\gamma }_i=1}{\overset{{\Pi }_i}{\Sigma }}(z_{s,k}^{\eta ,{\gamma }_i}-{\bar{z}}_{s,k}^{\eta }){(z_{s,k}^{\eta {,\gamma }_i}-{\bar{z}}_{s,k}^{\eta })}^{T,}$为机动目标状态采样粒子对应的预测观测值，模糊代价c_ζ,η可表示为对应的模糊隶属度观测值集合与机动目标状态采样粒子集合二维分配问题数学模型中的目标函数与约束条件可表示为：

目标函数$F\left(a\right)=\min \underset{\zeta =1}{\overset{{\Omega }_k}{\Sigma }}\underset{\eta =1}{\overset{{\Gamma }_k}{\Sigma }}{c}_{\zeta ,n}·a_{\zeta ,\eta }$

约束条件$\left(\begin{array}{cc}\underset{\eta =1}{\overset{{\Gamma }_k}{\Sigma }}{a}_{\zeta ,\eta }\le 1& \forall 1\le \zeta \le {\Omega }_k\\ \underset{\zeta =1}{\overset{{\Omega }_k}{\Sigma }}{a}_{\zeta ,\eta }\le 1& \forall 1\le \eta \le {\Gamma }_k\\ a_{\zeta ,\eta }=\mathrm{0,1}& 1\le \zeta \le {\Omega }_k,1\le \eta \le {\Gamma }_k\end{array}\right)$

其中，a_ζ,η为二进制决策变量，a_ζ,η＝1表示第ζ个观测值与第η个机动目标状态采样粒子相关联，a_ζ,η＝0表示其他，该二维分配问题的最优解可表示为每行、每列最多只有一个1的Ω_k×Γ_k阶关联决策矩阵

依据粒子群优化理论，定义粒子γ_pso的适应度值为其中，j_ζ,η为依据拍卖算法理论定义的竞拍价格，j_ζ,η＝λ₁g_ζ,η-λ₂v_ζ+ε，g_ζ,η为第ζ个观测值与第η个机动目标状态粒子集的关联价值，g_ζ,η＝max{c_ζ,η}-c_ζ,η，v_ζ为观测值的价值，λ₁为利益系数，λ₂为代价系数，ε为足够小的正数，设定粒子p_γpso为个体当前最优位置，p_g为群体中最佳粒子位置，粒子的位置和速度向量维数为Ω_k×Π_pso，Ω_k为观测集维数，Π_pso为粒子群的粒子数，第γ_pso个粒子的位置表示一组观测值与机动目标状态采样粒子的关联性判断结果，依据粒子群优化算法公式更新粒子的位置与速度，并计算每个粒子的适应度值，根据设定的迭代次数递归计算，则最优适应度值对应的竞拍价格即为全局最优值，相应粒子位置表示观测值集合与机动目标状态采样粒子集合的最优关联结果。对于机动目标是否新生与消失的判定问题，首先依据关联决策矩阵进行粗判断，若关联决策矩阵中某一行的元素值均近似为零，则认为该观测值可能来自于新生目标，或者来自于杂波；若关联决策矩阵中某一列的元素值均近似为零，则认为该目标可能消失。然后依据目标在观测场景中的存在概率进行精判断，对于 若P_e≥2/3，则认为时刻k新出现的观测值来自于新生目标t，此时在整个粒子集中增加该新生目标粒子集；对于若P_e≤1/3，则认为时刻k目 标t消失，此时在整个粒子集中删除该消失目标粒子集。

步骤5：在步骤4完成观测值集合与机动目标状态采样粒子集合关联判定的基础上，进行粒子权值计算，依据观测分辨单元内机动目标存在与不存在的似然比表达式 $L\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma }),L\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma })=\left(\begin{array}{cc}\frac{p\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma }=1)}{p\left(z_{s,k}\right|e_k^{t,\gamma }=0)}& e_k^{t,\gamma }=1\\ 1& e_k^{t,\gamma }=0\end{array}\right),$其中，为机动目标状态采样粒子存在的似然函数，为机动目标状态采样粒子不存在的似然函数。分三种情况进行处理，具体描述如下：

状态(a)：对于持续存在的机动目标采样粒子，依据预测粒子状态时判定出的目标模型状态及相应的目标状态转移密度函数，得持续存在的机动目标采样粒子γ_c的未归一化权值为${\tilde{w}}_k^{t,{\gamma }_c}=w_{k-1}^{t,\gamma }L\left(z_{s,k}\right|x_k^{t,{\gamma }_c},e_k^{t,{\gamma }_c}=1),$归一化权值为$w_k^{t,{\gamma }_c}=\frac{{\tilde{w}}_k^{t,{\gamma }_c}}{\underset{{\gamma }_i=1}{\overset{{\Pi }_c}{\Sigma }}{\tilde{w}}_k^{t,{\gamma }_i}}.$

状态(b)：对于新生的机动目标采样粒子，依据预测粒子状态时判定出的目标模型状态及相应的目标状态分布，得新生的机动目标采样粒子γ_b的非归一化权重为${\tilde{w}}_k^{t,{\gamma }_b}=\frac{L\left(z_{s,k}\right|x_k^{t,{\gamma }_b},{\tilde{e}}_k^{t,{\gamma }_b}=1)p\left(x_k^{t,{\gamma }_b}\right|{\tilde{e}}_k^{t,{\gamma }_b}=1,e_{k-1}^{t,\gamma }=0)}{{\Pi }_{b}q\left(x_k^{t,{\gamma }_b}\right|{\tilde{e}}_k^{t,{\gamma }_b}=1,e_{k-1}^{t,\gamma }=0,z_{s,k})},$归一化权值为$w_k^{t,{\gamma }_b}=\frac{{\tilde{w}}_k^{t,{\gamma }_b}}{\underset{{\gamma }_i=1}{\overset{{\Pi }_b}{\Sigma }}{\tilde{w}}_k^{t,{\gamma }_i}},$其中，为新生机动目标采样粒子先验分布，为新生机动目标采样粒子建议分布。

状态(c)：对于消失的机动目标采样粒子，粒子权值不存在。

在完成机动目标的粒子状态预测和粒子权重计算基础上，对持续存在的机动目标采样粒子与新生的机动目标采样粒子的粒子集进行混合。记机动目标新生概率为P_b，消失概率为P_d，k-1时刻机动目标存在概率为可得非归一化混合概率${\tilde{M}}_b=P_b[1-{\hat{P}}_{k-1}]\underset{{\gamma }_b=1}{\overset{{\Pi }_b}{\Sigma }}{\tilde{w}}_k^{t,{\gamma }_b},{\tilde{M}}_c=[1-P_d]{\hat{P}}_{k-1}\underset{{\gamma }_c=1}{\overset{{\Pi }_c}{\Pi }}{\tilde{w}}_k^{t,{\gamma }_c},$以及k时刻机动目标的存在概率更新持续存在的机动目标采样粒子与新生的机动目标采样粒子权重可得合并Π_b+Π_c个粒子，得k时刻机动目标混合粒子集其中i＝b,c，t＝0,1,…,Γ_k。步骤6：按照序贯重要性重采样理论，对步骤5获得的混合粒子集中的Π_b+Π_c个粒子进行重采样，得到逼近机动目标t状态后验概率分布的粒子集其中，k＝1,2,…，t＝0,1,…Γ_max，

步骤7：在时刻k，考虑粒子存在变量的影响，按照每个模型的模型概率进行粒子状态融合，得到机动目标t的局部状态后验估计值均方误差 $P_{k/k}^t=\underset{m_k^{t,\gamma }=1}{\overset{M}{\Sigma }}\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}{e}_k^{t,\gamma }{w}_k^{t,\gamma }({\hat{x}}_{k/k}^{m_k^{t,\gamma },t}-{\hat{x}}_{k/k}^t){({\hat{x}}_{k/k}^{m_k^{t,\gamma },t}-{\hat{x}}_{k/k}^t)}^T,$其中，对应于模型${\hat{x}}_{k/k}^{m_k^{t,\gamma },t}=\frac{\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}{e}_k^{t,\gamma }{x}_k^{t,\gamma }{w}_k^{t,\gamma }}{\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}{e}_k^{t,\gamma }},\gamma =\mathrm{1,2},...,\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }},$为预测模型概率，为对应模型的采样粒子数，${\mu }_k^{m_k^{t,\gamma }}=\frac{{\mu }_{k/k-1}^{m_k^{t,\gamma }}{L}_k^{m_k^{t,\gamma }}}{\underset{m_k^{t,\gamma }=1}{\overset{M}{\Sigma }}{\mu }_{k/k-1}^{m_k^{t,\gamma }}{L}_k^{m_k^{t,\gamma }}},$为模型似然，$L_k^{m_k^{t,\gamma }}=\frac{\underset{\gamma =1}{\overset{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}{\Sigma }}p\left(z_{s,k}\right|x_k^{t,\gamma },e_k^{t,\gamma }=1)}{\Pi {\mu }_{k/k-1}^{m_k^{t,\gamma }}}.$

依据目标存在的后验概率表达式P_e，得到时刻k机动目标数目Γ_k的估计值${\hat{\Gamma }}_k=\underset{t=1}{\overset{{\Gamma }_k}{\Sigma }}\left(\mathrm{Round}(P_e-\frac{1}{6})\right),$其中，$P_e\stackrel{\Delta }{=}P(e_k^{t,\gamma }=1|z_{s,k})=\frac{1}{\Pi }\underset{\gamma =1}{\overset{\Pi }{\Sigma }}{e}_k^{t,\gamma },$Round(x)表示求得距离x最近的整数。因此，在观测时间序列k＝1,2,…K内，对于传感器S_s，机动目标t的子航迹可表示为状态集合时刻k估计出的机动目标数目为

步骤8：在各时刻k，假设各传感器S_s之间已实现时空同步处理，各传感器信息处理系统之间的过程噪声、观测噪声互不相关，反复执行步骤3至步骤7，求得各传感器观测下机动目标t的S个局部状态估计值及其相应的估计均方误差阵 参照分布式融合结构，可得全局状态估计值全局状态估计均方误差观测场景中的目标数目估计值为各传感器局部估计值中的最大值，具体原理如图1所示。

为了更直观的表示本发明运用于可变数目机动目标跟踪时的效果，用本发明方法、经典多模粒子滤波算法(记方法一)与一种改进多模粒子滤波算法(记方法二)分别对图2中的机动目标进行跟踪。对于各机动目标，经过100次Monte Carlo仿真，图3给出了各时刻不同方法对典型目标的模型概率估计，图4给出了各时刻不同方法对典型目标的运动状态估计均方根误差，图5给出了各时刻不同方法对目标数目的估计。结果表明，在对目标运动模型的估计准确性、运动状态均方根误差估计与目标数目的估计方面，本发明方法的检测与跟踪性能均优于传统多模粒子滤波方法。这是由于本发明采用了粒子状态集与观测值的关联判定，当目标相距较近以及出现运动模型切换的时候，能够快速判定粒子集与观测值的关联度，使得多模型粒子滤波方法可以更有效地估计目标运动模型的概率，估计效果更接近目标运动的真实模型，能很好地适应目标的机动运动状态，进而做出较为准确的状态估计和目标数目估计；方法一和方法二由于没有考虑粒子集与观测值的关联判定，在低检测率场景中，当目标运动状态相近或运动模型切换的时候会导致目标运动模型的错误估计，容易出现误跟踪及目标数目估计不准确。当目标新出现或消失时，本发明所采用的软判定准则不能及时判定出目标数目变化，一定程度上影响了运动状态估计精度，但变化间隔越短所受影响越小；而方法一和方法二的硬性判定虽然可以相对快速地检测到目标数目变化，但是当受较强杂波影响时会影响目标数目的准确估计，进而导致目标状态的错误估计。因此，在低信噪比条件下，本发明方法能有效处理目标相距较近或者目标数目发生变化的机动目标检测与跟踪。