CN107220993B - 基于Mean shift聚类的粒子滤波机动目标追踪方法 - Google Patents

Abstract

本发明涉及一种基于Mean shift聚类的粒子滤波机动目标追踪算法，该算法首先利用Mean shift的特征空间方法自适应的计算粒子集的聚类，无需初始化聚类中心与聚类个数，便能自适应的准确确定粒子集中的聚类，通过聚类精确表示目标。其次在进行目标概率密度估计与相似性度量时，使用核函数描述像素点的位置信息，综合像素点的位置信息与空间信息进行目标概率估计，精确描述目标模板，增加目标追踪精度。在粒子滤波追踪结果确定时，对当前追踪结果概率模型赋予一定权值，对目标模板进行更新，可以增强算法的动态性，提高追踪过程的抗干扰能力。

Description

基于Mean shift聚类的粒子滤波机动目标追踪方法

技术领域

本发明涉及一种空间机动目标的追踪方法，特别是一种基于Mean shift聚类的粒子滤波机动目标追踪算法。

背景技术

文献“基于遗传算法的改进粒子滤波算法，上海交通大学学报，2011，Vol45(10)，p1526-1530”提出了一种基于遗传算法的改进粒子滤波目标追踪算法。该方法针对粒子数量和质量造成的粒子滤波中粒子的退化问题，采用遗传算法进行初始粒子优化，改善初始粒子质量。通过顺序重要性采样(SIS)获得一组大量粒子，再通过遗传算法优化出一组数量较少的粒子作为粒子滤波的初始粒子。该追踪算法既能改善初始粒子质量，又能保证其随机性和统计性，减少粒子退化的可能性。文献所述追踪算法将遗传算法引入粒子重采样中，用来增加粒子多样性；但优选方法单一，粒子多样性小，新生粒子精度低，会造成计算量急剧增加，追踪速率低，不满足工业实时性要求。该算法在进行目标模板描述时，没有综合考虑颜色信息与空间信息对模板相似性的影响，影响算法追踪结果的准确性，且目标模板仅依赖于初始帧目标的选择，当目标发生形变等情况时，目标模板无相应改变，因而该方法不具有动态性，抗干扰能力弱。

发明内容

要解决的技术问题

为了克服现有的基于遗传算法的改进粒子滤波目标追踪算法的准确性、实时性与动态性的不足，本发明提供一种基于Mean shift聚类的粒子滤波机动目标追踪算法。

技术方案

一种基于Mean shift聚类的粒子滤波机动目标追踪算法，其特征在于步骤如下：

步骤1：定义初始时刻的粒子集为S₁，粒子集S₁包含粒子个数为N；每个粒子模型为：

式中，X^j为粒子集中第j个粒子，为粒子中心位置坐标，分别为粒子在x、y方向的运动速度，分别为粒子所在区域的宽和高，a是尺度因子；

步骤2：以目标模板的中心点为基准点，采用核函数计算各点的权重，则对目标模板进行概率密度估计为：

式中，q_u表示目标模板概率密度直方图，h为核函数带宽，δ[b(x_i)-u]为直方图索引函数，u为直方图中任一直方块的编号，M为目标模板内像素点的个数，x_i为模板内任一像素点，x为模板中心像素点，C为归一化系数，为核函数。

步骤3：判断当前时刻是否为初始时刻，若是初始时刻，转入步骤4；若不是初始时刻，则对粒子集进行重采样，已知t-1时刻的传播粒子集为计算粒子集中各粒子的归一化权值累积量：

式中，为t-1时刻粒子集中第j个粒子的归一化权值累积量，为t-1时刻第i个粒子的归一化权值，为t-1时刻经Mean shift聚类后粒子集包含的粒子数；

生成N个在[0,1]区间均匀分布的随机数对于每一随机数μ^j，寻找归一化权值累计量大于或等于λ^j的最小标号m，即λ_m-1＜μ^j＜λ_m；当μⁱ落在区间[λ_m-1,λ_m]时， t-1时刻传播粒子集的第m个粒子被复制一次，N次复制后，获得更新的包含 N个粒子的样本集合S_t；

步骤4：利用Mean shift对粒子集S_t进行聚类，对每一粒子进行多次迭代，计算各粒子对应的聚类中心；对于粒子集S_t中任一粒子X^j，其第k+1次迭代时Mean shift向量为：

式中，h为带宽，X^j,k为粒子集S_t中第j个粒子的第k次迭代获得的粒子，P是以粒子X^j,k位置为圆心、以h为直径的圆内包含的除该粒子外其他粒子个数；为核函数；

Mean shift向量表示粒子每次迭代前后粒子位置的变化，则每次Mean shift迭代后粒子位置计算公式为：

式中，X^j,k为粒子集S_t中第j个粒子的第k次迭代获得的粒子，X^j,k+1为粒子集S_t中第j个粒子的第k+1次迭代获得的粒子；

每迭代一次，迭代结果更接近聚类中心，迭代次数越多，每次迭代位置变化的距离会越少；用||X^j,k+1-X^j,k||表示粒子X^j,k+1位置与粒子X^j,k的位置之间的距离，当 ||X^j,k+1-X^j,k||≤ε，ε为10像素时，停止迭代，此时认为X^j,k+1为粒子X^j对应的聚类中心，迭代计算出粒子集S_t中所有粒子对应的聚类中心，迭代收敛至同一位置的粒子属于同一聚类中心；所有聚类中心点组成新的经Mean shift聚类的粒子集此时粒子集总数从N变为

步骤5：根据粒子传播方程确定粒子集中各个粒子下一时刻将传播的位置，获得对应的传播粒子集粒子传播状态方程为：

式中，为粒子集中任一粒子，为粒子对应的传播粒子，且为粒子集中的粒子，w_t为高斯噪声，A为粒子传播矩阵；

据此可获得粒子集对应的传播粒子集粒子集包含粒子数仍为

以粒子为中心的样本模板概率密度与目标模板概率密度的Bhattacharyya系数为：

式中，q_u，t为t时刻目标模板的概率密度直方图，表示t时刻粒子集中第j个粒子，表示以为中心的样本模板的概率密度直方图，u为色彩直方图直方块的编号，m为色彩直方图u的个数；

将以粒子为中心的样本模板概率密度与目标模板概率密度的Bhattacharyya系数进行归一化，计算出粒子的归一化权值，并将其作为对应的粒子的归一化权值，则粒子的归一化权值为：

式中，表示t时刻粒子集中第j个粒子的归一化权值，为t时刻经过Meanshift聚类后粒子集包含的粒子个数；

则目标跟踪输出为：

式中，表示t时刻粒子集中第j个粒子的位置，表示粒子对应的归一化权值，为t时刻经过Mean shift聚类后粒子集包含粒子个数，E[S_t]为t时刻目标跟踪输出结果；

步骤6：更新追踪结果区域的概率密度直方图q_u,t：

q_u,t+1＝0.8×q_u,t+0.2×p_u,t

式中，q_u,t为t时刻目标模板的概率密度直方图，q_u,t+1为t+1时刻目标模板的概率密度直方图，p_u,t为t时刻跟踪输出位置的概率密度直方图；

步骤7：循环步骤3～6，进行下一时刻的目标追踪。

有益效果

本发明提出的一种基于Mean shift聚类的粒子滤波机动目标追踪算法，本算法在重采样之后，对粒子集进行Mean shift聚类分析，算法无需对聚类中心及聚类数量进行初始化，能自适应的准确计算特征空间中的聚类，由此产生的预测粒子可较好地服从状态的后验概率分布，从而有效地减少描述后验概率密度函数所需的粒子数，既简化了数据复杂度，提高目标追踪效率，又保证了追踪准确性。加入了空间信息描述目标模板概率分布，使目标追踪的准确性得到提高。当每一时刻获得跟踪输出位置后，对目标模板进行更新，可在目标发生形变，旋转等情况下保证算法的追踪效果。

具体实施方式

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

该算法首先利用Mean shift的特征空间方法自适应的计算粒子集的聚类，无需初始化聚类中心与聚类个数，便能自适应的准确确定粒子集中的聚类，通过聚类精确表示目标。其次在进行目标概率密度估计与相似性度量时，使用核函数描述像素点的位置信息，综合像素点的位置信息与空间信息进行目标概率估计，精确描述目标模板，增加目标追踪精度。在粒子滤波追踪结果确定时，对当前追踪结果概率模型赋予一定权值，对目标模板进行更新，可以增强算法的动态性，提高追踪过程的抗干扰能力。

1、粒子样本集合初始化。定义初始时刻的粒子集为S₁，粒子集S₁包含粒子个数为N，粒子模型为：

式中，X^j为粒子集中第j个粒子，为粒子中心位置坐标，分别为粒子在x、y方向的运动速度，分别为粒子所在区域的宽和高，a是尺度因子。

2、目标模板概率密度描述。

加入像素点的空间位置信息，以目标模板的中心点为基准点，依据各像素点与基准点的距离，引进核函数为不同像素点赋予不同权重，采用以下核函数计算各像素点的权重：

式中，x_mid为目标模板的x轴中心坐标，y_mid为目标模板的y轴中心坐标。x_i,y_i分别点x 的横、纵坐标，(2*x_mid)²+(2*y_mid)²是核函数半径的平方。当点x在目标模板中心时， K(x)的值为1，当点x在模板边界时，K(x)的值为3/4，距基准点越远的像素点对应的权值越小。

则对目标模板进行概率密度估计为：

式中，q_u表示目标模板概率密度直方图，h为核函数带宽，δ[b(x_i)-u]为直方图索引函数，u为直方图中任一直方块的编号，M为给定空间内像素点的个数，x_i为模板内任一像素点，C为归一化系数，为核函数。

3、对粒子集进行重采样。

判断当前时刻是否为初始时刻，若是初始时刻，转入步骤4，若不是初始时刻，则对粒子集进行重采样，已知t-1时刻的传播粒子集为计算粒子集现中各粒子的归一化权值累积量：

式中，为t-1时刻粒子集中第j个粒子的归一化权值累积量，为t-1时刻粒子集第j个粒子的归一化权值，为粒子集包含的粒子数。

生成N个在[0，1]区间均匀分布的随机数构成数组r。对于每一随机数μ^j，分别用二分查找法搜索归一化累计概率数组，令μ^j满足：

λ_m-1＜μ^j＜λ_m (5)

式中，λ_m-1为粒子集中第m-1个粒子的累计归一化权值，λ_m为粒子集中第m个粒子的累计归一化权值，则m为随机数μ^j对应的归一化权值累计量大于或等于μ^j的最小标号，即当μ^j落在区间[λ_m-1,λ_m]时，被复制一次，N次复制后，获得更新的粒子集S_t，粒子集S_t包含N个粒子。

4、利用Mean shift对粒子集进行聚类。

利用Mean shift对粒子集S_t进行聚类，对每一粒子进行多次迭代，计算各粒子对应的聚类中心。对于任一粒子X^j，其第k+1次迭代时Mean shift向量为：

式中，h为带宽，X^j,k为粒子集S_t中第j个粒子的第k次迭代获得的粒子，P是以粒子X^j,k位置为圆心、以h为半径的圆内包含的除该粒子外其他粒子个数。Xⁱ代表以粒子X^j,k位置为圆心、以h为直径的圆内包含的除该粒子其他任意粒子。

Mean shift向量表示粒子每次迭代前后粒子位置的变化，则每次Mean shift迭代后粒子位置计算公式为：

式中，X^j,k为粒子集S_t中第j个粒子的第k次迭代获得的粒子，X^j,k+1为粒子集S_t中第j个粒子的第k+1次迭代获得的粒子。

每迭代一次，迭代结果更接近聚类中心，迭代次数越多，每次迭代位置变化的距离会越少。用||X^j,k+1-X^j,k||表示X^j,k+1与X^j,k之间的距离，当||X^j,k+1-X^j,k||≤ε，ε为10 像素时，停止迭代，此时认为X^j,k+1为粒子X^j对应的聚类中心，迭代计算出粒子集S_t中所有粒子对应的聚类中心，迭代收敛至同一位置的粒子属于同一聚类中心。所有聚类中心点组成新的经Mean shift聚类的粒子集此时粒子集总数从N变为

5、粒子传播与目标位置预估。

根据粒子传播方程确定粒子集中各个粒子下一时刻将传播的位置，获得对应的传播粒子集粒子传播状态方程为：

式中，为粒子集中任一粒子，为粒子对应的传播粒子，且粒子属于传播粒子集w_t为高斯噪声，A为粒子传播矩阵，具体表示如下：

式中，为粒子中心位置坐标，分别为粒子在x、y方向的运动速度，分别为粒子所在区域的宽和高，是尺度因子。

式中，为粒子中心位置坐标，分别为粒子在x、y方向的运动速度，分别为粒子所在区域的宽和高，是尺度因子。

式中，A为粒子传播矩阵，分别为粒子所在区域的宽和高，ΔT表示粒子传播的时间间隔。

据此可获得粒子集对应的传播粒子集粒子集包含粒子数仍为

计算以粒子为中心的模板概率密度估计：

式中，p_u表示以粒子为中心的样本模板的概率密度直方图，h为核函数带宽，δ[b(x_i)-u]为直方图索引函数，u为直方图中任一直方块的编号，M为样本模板内像素点的个数，x_i为模板内任一像素点，x为模板中心处的像素点，即粒子位置处的像素点，C为归一化系数。

以粒子为中心的样本模板概率密度与目标模板概率密度的Bhattacharyya系数为：

式中，q_u,t为t时刻目标模板的概率密度直方图，表示t时刻粒子集中第j个粒子，表示以为中心的样本模板的概率密度直方图，u为色彩直方图直方块的编号，m为色彩直方图u的个数。

将以粒子为中心的模板概率密度与目标模板概率密度的Bhattacharyya系数，进行归一化，计算出粒子的归一化权值，并将其作为对应的粒子的归一化权值，则粒子的归一化权值为：

式中，表示t时刻粒子集中第j个粒子的归一化权值，为t时刻经过Meanshift 聚类后粒子集包含的粒子个数。

则目标跟踪输出为：

式中，表示t时刻粒子集中第j个粒子的位置，表示粒子对应的归一化权值，为经过Mean shift聚类后粒子集包含粒子个数，E[S_t]为t时刻目标跟踪输出结果。

6、模型更新。

计算追踪结果区域的概率密度直方图q_u,t：

式中，q_u,t为t时刻追踪结果区域的概率密度直方图，h为核函数带宽，δ[b(x_i)-u]为直方图索引函数，u为直方图中任一直方块的编号，M为追踪模板内像素点的个数，x_i为追踪结果区域内任一像素点，x为追踪结果区域的中心像素点，C为归一化系数，为核函数。

对下一时刻的追踪目标模板进行更新：

q_u,t+1＝0.8×q_u,t+0.2×p_u,t (17)

式中，q_u,t为t时刻目标模板的概率密度直方图，q_u,t+1为t+1时刻目标模板的概率密度直方图，p_u,t为t时刻跟踪输出位置的概率密度直方图。

循环步骤3～步骤6，进行下一时刻的目标追踪。

Claims (1)

1.一种基于Mean shift聚类的粒子滤波机动目标追踪方法，其特征在于步骤如下：

步骤1：定义初始时刻的粒子集为S₁，粒子集S₁包含粒子个数为N；每个粒子模型为：

式中，X^j为粒子集中第j个粒子，为粒子中心位置坐标，分别为粒子在x、y方向的运动速度，分别为粒子所在区域的宽和高，a是尺度因子；

步骤2：以目标模板的中心点为基准点，采用核函数计算各点的权重，则对目标模板进行概率密度估计为：

式中，q_u表示目标模板概率密度直方图，h为核函数带宽，δ[b(x_i)-u]为直方图索引函数，u为直方图中任一直方块的编号，M为目标模板内像素点的个数，x_i为模板内任一像素点，x为模板中心像素点，C为归一化系数，为核函数；

步骤3：判断当前时刻是否为初始时刻，若是初始时刻，转入步骤4；若不是初始时刻，则对粒子集进行重采样，已知t-1时刻的传播粒子集为计算粒子集中各粒子的归一化权值累积量：

式中，为t-1时刻粒子集中第j个粒子的归一化权值累积量，为t-1时刻第i个粒子的归一化权值，为t-1时刻经Mean shift聚类后粒子集包含的粒子数；

生成N个在[0,1]区间均匀分布的随机数对于每一随机数μ^j，寻找归一化权值累计量大于或等于λ^j的最小标号m，即λ_m-1＜μ^j＜λ_m；当μ^j落在区间[λ_m-1,λ_m]时，t-1时刻传播粒子集的第m个粒子被复制一次，N次复制后，获得更新的包含N个粒子的样本集合S_t；

步骤4：利用Mean shift对粒子集S_t进行聚类，对每一粒子进行多次迭代，计算各粒子对应的聚类中心；对于粒子集S_t中任一粒子X^j，其第k+1次迭代时Mean shift向量为：

式中，h为带宽，X^j,k为粒子集S_t中第j个粒子的第k次迭代获得的粒子，P是以粒子X^j,k位置为圆心、以h为直径的圆内包含的除该粒子外其他粒子个数；为核函数；Xⁱ代表以粒子X^j,k位置为圆心、以h为直径的圆内包含的除该粒子的 其他任意粒子；

Mean shift向量表示粒子每次迭代前后粒子位置的变化，则每次Mean shift迭代后粒子位置计算公式为：

式中，X^j,k为粒子集S_t中第j个粒子的第k次迭代获得的粒子，X^j,k+1为粒子集S_t中第j个粒子的第k+1次迭代获得的粒子；

每迭代一次，迭代结果更接近聚类中心，迭代次数越多，每次迭代位置变化的距离会越少；用||X^j,k+1-X^j,k||表示粒子X^j,k+1位置与粒子X^j,k的位置之间的距离，当||X^j,k+1-X^j,k||≤ε，ε为10像素时，停止迭代，此时认为X^j,k+1为粒子X^j对应的聚类中心，迭代计算出粒子集S_t中所有粒子对应的聚类中心，迭代收敛至同一位置的粒子属于同一聚类中心；所有聚类中心点组成新的经Mean shift聚类的粒子集此时粒子集总数从N变为

步骤5：根据粒子传播方程确定粒子集中各个粒子下一时刻将传播的位置，获得对应的传播粒子集粒子传播状态方程为：

式中，为粒子集中任一粒子，为粒子对应的传播粒子，且为粒子集中的粒子，w_t为高斯噪声，A为粒子传播矩阵；

据此可获得粒子集对应的传播粒子集粒子集包含粒子数仍为

以粒子为中心的样本模板概率密度与目标模板概率密度的Bhattacharyya系数为：

式中，q_u,t为t时刻目标模板的概率密度直方图，表示t时刻粒子集中第j个粒子，表示以为中心的样本模板的概率密度直方图，u为直方图中任一直方块的编号，m为色彩直方图u的个数；

将以粒子为中心的样本模板概率密度与目标模板概率密度的Bhattacharyya系数进行归一化，计算出粒子的归一化权值，并将其作为对应的粒子的归一化权值，则粒子的归一化权值为：

式中，表示t时刻粒子集中第j个粒子的归一化权值，为t时刻经过Meanshift聚类后粒子集包含的粒子个数；

则目标跟踪输出为：

式中，表示t时刻粒子集中第j个粒子的位置，表示粒子对应的归一化权值，为t时刻经过Mean shift聚类后粒子集包含粒子个数，E[S_t]为t时刻目标跟踪输出结果；

步骤6：计算追踪结果区域的概率密度直方图q_u,t：

式中，q_u,t为t时刻追踪结果区域的概率密度直方图，h为核函数带宽，δ[b(x_i)-u]为直方图索引函数，u为直方图中任一直方块的编号，M为追踪模板内像素点的个数，x_i为追踪结果区域内任一像素点，x为追踪结果区域的中心像素点，C为归一化系数，为核函数；

对下一时刻的追踪目标模板q_u,t+1进行更新：

q_u,t+1＝0.8×q_u,t+0.2×p_u,t

式中，q_u,t为t时刻目标模板的概率密度直方图，q_u,t+1为t+1时刻目标模板的概率密度直方图，p_u,t为t时刻跟踪输出位置的概率密度直方图；

步骤7：循环步骤3～6，进行下一时刻的目标追踪。