CN104215960A - 基于改进粒子滤波的目标跟踪方法 - Google Patents

Abstract

本发明属于雷达目标跟踪技术领域，公开了基于改进粒子滤波的目标跟踪方法。该基于改进粒子滤波的目标跟踪方法包括以下步骤：对回波数据进行脉压处理，得到K帧观测数据并确定相应的恒虚警检测门限，根据检测门限确定第k帧的目标影响区域C_k；当k＝1时，产生N个初始粒子，抽取N_c个继续粒子进行重采样，得出N_c个重采样粒子以及对应的权值，计算出目标状态值x₁；当k>1时，利用第k-1帧的N_c个重采样粒子得到N_c个预测粒子并抽取N_b个出生粒子，得到第k帧的N个粒子；计算第k帧粒子的权值并进行归一化处理；抽取第k帧的继续粒子并对继续粒子进行重采样，计算目标状态值x_k；根据每一帧的结果得到检测航迹。本发明与现有技术相比，计算量小，粒子多样性增加，可用于雷达微弱目标检测与跟踪。

Description

基于改进粒子滤波的目标跟踪方法

技术领域

本发明属于雷达目标跟踪技术领域，特别涉及基于改进粒子滤波的目标跟踪方法。

背景技术

目前由于检测前跟踪算法(Track Before Detect,TBD)算法具有优良的弱小目标检测能力，逐步被引入到雷达目标检测中。在TBD算法中，粒子滤波算法结合蒙特卡罗积分与贝叶斯估计，利用一系列随机抽取的样本以及样本的权重来表示状态的后验概率分布，相比于传统的卡尔曼滤波等算法具有不需要满足系统为线性，噪声为高斯分布，后验概率也是高斯型等限制条件的优点，在雷达、声纳、红外、导航制导、计算机视觉、金融统计与目标跟踪等领域得到了广泛应用。

虽然粒子滤波具有良好的性能，对研究微弱目标检测具有很大贡献，然而粒子滤波算法仍未成熟，它还存在着以下的一些问题：1)计算量太大。这是粒子滤波用于工程实践的最大障碍。主要有两个原因：一是当选择的重要性函数偏离真实的状态后验概率较多时，要获得较高的滤波精度，必须在状态空间分布大量的随机样本进行搜索；二是重采样过程增加了算法的复杂度。Jayesh和Peter提出了高斯粒子滤波算法，该算法与标准粒子滤波相比计算量大大减少，复杂度降低，但是高斯滤波在后验分布不能用高斯分布近似的非线性动态空间模型或非线性系统非高斯噪声模型时，滤波性能不佳。2)应用于复杂问题时粒子滤波会变得不稳定。统计表明，即使粒子滤波产生了适当的权重粒子，基于这些粒子的随机变量预测也可能溢出而使得预测不可靠。这是因为，粒子滤波的性能会由于蒙特卡洛随机采样的随机特性而相差非常大，这在滤波过程中尤其有害，因为误差可能在滤波过程中成指数关系进行积累，甚至使粒子滤波器失效。3)粒子容易丧失多样性。在重采样过程中，权值较大的粒子被多次选中，而权值较小的粒子被逐渐抛弃，最后所有粒子都崩溃到同一个点上，严重丧失了粒子多样性，粒子滤波的估计性能降低，在系统噪声方差较小的这种情况下尤为严重。当前比较常见的马尔科夫蒙特卡罗(MCMC)算法与标准粒子滤波算法相比粒子多样性增加，但是该算法的计算量也明显增加。

发明内容

本发明针对现有标准粒子滤波算法计算量太大，蒙特卡洛随机采样产生的粒子不均匀以及重采样造成的粒子多样性丧失的问题，提出了基于改进粒子滤波的目标跟踪方法，与标准粒子滤波方法相比，本发明的改进粒子滤波方法计算量较少，粒子多样性增加。

为实现上述技术目的，本发明采用如下技术方案予以实现。

基于改进粒子滤波的目标跟踪方法包括以下步骤：

步骤1，当目标在一个平面内运动时，利用雷达接收目标的回波数据；对目标的回波数据进行距离向脉冲压缩处理，得到第1帧观测数据至第K帧观测数据，K为自然数；

步骤2，针对每帧观测数据，确定恒虚警检测门限；

步骤3，将每帧观测数据中低于恒虚警检测门限的数据丢弃，保留每帧观测数据中大于或等于恒虚警检测门限的数据；将第k帧观测数据中保留数据对应的区域记为第k帧观测数据的目标影响区域C_k，k取1至K；

步骤4，当k＝1时，在第k帧观测数据的目标影响区域C_k用Halton拟随机序列均匀地产生N个初始粒子，N为设定的自然数；得出第1帧观测数据对应的每个初始粒子的权值；对第1帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第1帧观测数据对应的每个初始粒子的归一化权值；在第1帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个初始粒子作为第1帧观测数据对应的继续粒子，N_c为自然数且N_c<N；对第1帧观测数据对应的每个继续粒子进行重采样，得出第1帧观测数据对应的N_c个重采样粒子、以及第1帧观测数据对应的每个重采样粒子的权值；第1帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第1帧观测数据对应的每个重采样粒子的权值为计算出第1帧观测数据对应的目标状态值x₁：

$x_1=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_1^{\left(m\right)}{x}_1^{\left(m\right)};$

步骤5，当k>1时，将第k-1帧观测数据对应的第m个重采样粒子代入重要性函数q(x_k|x_k-1,z_k)中进行采样，得到第k帧观测数据对应的第m个预测粒子将m从1至N_c进行遍历，得出第k帧观测数据对应的N_c个预测粒子；按均匀分布密度函数抽取N_b个出生粒子，N_b＝N-N_c；

步骤6，得出第k帧观测数据对应的N个初始粒子，第k帧观测数据对应的N个初始粒子包括步骤5中第k帧观测数据对应的N_c个预测粒子和N_b个出生粒子；得出第k帧观测数据对应的每个初始粒子的权值；

步骤7，对第k帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第k帧观测数据对应的每个初始粒子的归一化权值；

步骤8，在第k帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第k帧观测数据对应的继续粒子；

步骤9，对第k帧观测数据对应的每个继续粒子进行重采样，得出第k帧观测数据对应的N_c个重采样粒子、以及第k帧观测数据对应的每个重采样粒子的权值；第k帧观测数据对应的第m个重采样粒子表示为第k帧观测数据对应的每个重采样粒子的权值为计算出第k帧观测数据对应的目标状态值x_k：

$x_k=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_k^{\left(m\right)}{x}_k^{\left(m\right)}.$

本发明的特点和进一步改进在于：

在步骤1中，以目标的运动平面为基准建立二维平面直角坐标系，二维平面直角坐标系的横轴为X轴，纵轴为Y轴；

第k帧观测数据表示为z_k， $z_k=\left\{z_k^{(i,j)}\right|1\&le;i\&le;M,1\&le;j\&le;N,1\&le;k\&le;K\},$ 表示第k帧观测数据中位于第i行第j列的像素单元对应的观测数据；M表示每帧观测数据对应的像素单元的行数，N表示每帧观测数据对应的像素单元的列数。

所述步骤4的具体子步骤为：

(4.1)当k＝1时，在第k帧观测数据的目标影响区域C_k用Halton拟随机序列均匀地产生N个初始粒子的过程为：将十进制数dec1分别转化为2进制数和3进制数，dec1为自然数且dec1的取值为1至N；分别对2进制数和3进制数求反序，得到2进制反序后数和3进制反序后数；在2进制反序后数依次前置一个小数点和0，得到2进制小数，同时，在3进制反序后数依次前置一个小数点和0，得到3进制小数；将2进制小数转化为十进制数dec2，将3进制小数转化为十进制数dec3；则产生的第dec1个初始粒子对应X方向上的第round(l_X,1·dec2))个位置，同时对应Y方向上的第round(l_Y,1·dec3))个位置，round(·)表示取四舍五入，l_X,1表示第1帧观测数据的目标影响区域C₁中X方向的长度，l_Y,1表示第1帧观测数据的目标影响区域C₁中Y方向的长度；

(4.2)得出第1帧观测数据对应的每个初始粒子的权值，第1帧观测数据对应的第n个初始粒子的权值为：

其中，表示第1帧观测数据对应的第n个初始粒子值，表示第1帧观测数据的目标影响区域C₁中像素单元行数的集合，表示第1帧观测数据的目标影响区域C₁中像素单元列数的集合；为：

$l\left(z_1^{(i,j)}\right|{\hat{x}}_1^{\left(n\right)})=\exp [-\frac{h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)(h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)-2z_1^{(i,j)})}{{2\mathrm{\&sigma;}}_n^2}]$

其中，是观测噪声方差，的含义，表示坐标为(x₁,y₁)，强度为I₁的目标对位于第i行第j列的像素单元的强度贡献，其表达式为：

$h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)=\frac{\mathrm{\&Delta;x\&Delta;y}{I}_1}{2\mathrm{\&pi;}{\mathrm{\&Gamma;}}^2}\exp [-\frac{{(\mathrm{i\&Delta;x}-x_1)}^2+{(\mathrm{j\&Delta;y}-y_1)}^2}{{2\mathrm{\&Gamma;}}^2}]$

其中，Δx表示每帧观测数据对应的X方向分辨率，Δy表示每帧观测数据对应的Y方向分辨率，Γ为已知的雷达模糊斑点数量；

(4.3)对第1帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第1帧观测数据对应的每个初始粒子的归一化权值；第1帧观测数据对应的第n个初始粒子的归一化权值为：

$w_1^{\left(n\right)}=\frac{{\tilde{w}}_1^{\left(n\right)}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\tilde{w}}_1^{\left(n\right)}}$

其中，n取1至N；

(4.4)在第1帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第1帧观测数据对应的继续粒子；

(4.5)第1帧观测数据对应的每个继续粒子进行重采样，得出第1帧观测数据对应的N_c个重采样粒子、以及第1帧观测数据对应的每个重采样粒子的权值；第1帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第1帧观测数据对应的每个重采样粒子的权值为计算出第1帧观测数据对应的目标状态值x₁：

$x_1=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_1^{\left(m\right)}{x}_1^{\left(m\right)}.$

在步骤6中，第k帧观测数据对应的第n个初始粒子的权值为：

其中，表示第k帧观测数据对应的第n个初始粒子值，表示第k帧观测数据的目标影响区域C_k中像素单元行数的集合，表示第k帧观测数据的目标影响区域C_k中像素单元列数的集合；为为：

$l\left(z_k^{(i,j)}\right|{\hat{x}}_k^{\left(n\right)})=\exp [-\frac{h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)(h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)-2z_k^{(i,j)})}{{2\mathrm{\&sigma;}}_n^2}]$

其中，是观测噪声方差，表示坐标为(x_k,y_k)，强度为I_k的目标对位于第i行第j列的像素单元的强度贡献，其表达式为：

$h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)=\frac{\mathrm{\&Delta;x\&Delta;y}{I}_k}{2\mathrm{\&pi;}{\mathrm{\&Gamma;}}^2}\exp [-\frac{{(\mathrm{i\&Delta;x}-x_k)}^2+{(\mathrm{j\&Delta;y}-y_k)}^2}{{2\mathrm{\&Gamma;}}^2}]$

其中，Δx表示每帧观测数据对应的X方向分辨率，Δy表示每帧观测数据对应的Y方向分辨率；Γ为已知的雷达模糊斑点数量。

在步骤7中，第k帧观测数据对应的第n个初始粒子的归一化权值为：

$w_k^{\left(n\right)}=\frac{{\tilde{w}}_k^{\left(n\right)}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\tilde{w}}_k^{\left(n\right)}}$

其中，n取1至N，为第k帧观测数据对应的第n个初始粒子的权值。

在步骤9之后，根据每帧观测数据对应的目标状态值，得出目标航迹。

本发明的有益效果为：1)本发明利用恒虚警设置低门限对回波数据预处理，只保留超过门限的数据点，使粒子尽量分布在受目标影响的较小区域，减小搜索范围，从而减少了计算量。2)本发明利用Halton拟随机序列在目标影响区域均匀的产生粒子，避免粒子产生群聚现象，使粒子分布的更加均匀，精度更高。3)本发明还将粒子划分为继续粒子和自适应出生粒子，主要思想就是保留权值较大的继续粒子进行动态更新，淘汰权值较小的粒子，同时根据粒子权值的大小分布情况，自适应地产生出生粒子，从而增加了粒子多样性。

附图说明

图1为本发明的基于改进粒子滤波的目标跟踪方法的流程图；

图2为环境信噪比为3dB时两种方法得出的检测结果比较示意图；

图3为环境信噪比为3dB时两种方法得出的距离检测均方误差比较示意图；

图4为环境信噪比为5dB时两种方法得出的检测结果比较示意图；

图5为环境信噪比为3dB时两种方法得出的距离检测均方误差比较示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的基于改进粒子滤波的目标跟踪方法的流程图。该基于改进粒子滤波的目标跟踪方法包括以下步骤：

步骤1，当目标在一个平面内运动时，利用雷达接收目标的回波数据。以目标的运动平面为基准建立二维平面直角坐标系，二维平面直角坐标系的横轴为X轴，纵轴为Y轴。然后，对目标的回波数据进行距离向脉冲压缩处理，得到第1帧观测数据至第K帧观测数据。其中，第k帧观测数据表示为z_k， $z_k=\left\{z_k^{(i,j)}\right|1\&le;i\&le;M,1\&le;j\&le;N,1\&le;k\&le;K\},$ 表示第k帧观测数据中位于第i行第j列的像素单元对应的观测数据；M表示每帧观测数据对应的像素单元的行数，N表示每帧观测数据对应的像素单元的列数。

步骤2，针对每帧观测数据，确定恒虚警检测门限；

其具体子步骤为：

利用基于自动删除算法的最小选择恒虚警率处理方法(参见“基于自动删除算法的最小选择恒虚警检测器”，许江湖张明敏王平波，2006,12)，在检测单元两侧的前、后沿滑窗中同时利用二元假设检验方法自动删除可能存在的干扰目标样本，然后分别求两个滑窗中剩余样本的平均值，取最小者作为噪声背景功率水平的估计，从而确定恒虚警检测门限(自适应检测门限)。

步骤3，对于每帧观测数据作低门限恒虚警数据预处理，将每帧观测数据中低于恒虚警检测门限的数据丢弃，保留每帧观测数据中超过(大于或等于)恒虚警检测门限的数据。将第k帧观测数据中保留数据对应的区域记为第k帧观测数据的目标影响区域C_k。

步骤4，当k＝1时，在第k帧观测数据的目标影响区域C_k用Halton拟随机序列均匀地产生N个初始粒子，N为设定的自然数；得出第1帧观测数据对应的每个初始粒子的权值；对第1帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第1帧观测数据对应的每个初始粒子的归一化权值；在第1帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第1帧观测数据对应的继续粒子，N_c为自然数且N_c<N；对第1帧观测数据对应的每个继续粒子进行重采样，得出第1帧观测数据对应的N_c个重采样粒子、以及第1帧观测数据对应的每个重采样粒子的权值；第1帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第1帧观测数据对应的每个重采样粒子的权值为计算第1帧观测数据对应的目标状态值x₁：

$x_1=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_1^{\left(m\right)}{x}_1^{\left(m\right)}.$

其具体子步骤为：

(4.1)当k＝1时，在第k帧观测数据的目标影响区域C_k用Halton拟随机序列均匀地产生N个初始粒子，N为设定的自然数。具体地，Halton拟随机序列的产生机制为：将十进制数dec1(dec1为自然数，其取值为1至N)分别转化为2进制数和3进制数(如十进制数6转化为2机制数110)，分别对2进制数和3进制数求反序，得到2进制反序后数和3进制反序后数，例如，2机制数110求反序后变为2进制反序后数011。在2进制反序后数依次前置一个小数点和0，得到2进制小数，例如，2进制反序后数011在依次前置一个小数点和0，变为2进制小数0.011。同时，在3进制反序后数依次前置一个小数点和0，得到3进制小数。将2进制小数转化为十进制数dec2(例如，2进制小数0.011转化为十进制数0.375)，将3进制小数转化为十进制数dec3。由于dec1的取值为1至N，则dec2为对应的N个取值在[0,1]中的值，同样，dec3为对应的N个取值在[0,1]中的值。此时，第dec1个初始粒子对应X方向上的第round(l_X,1·dec2))个位置，同时对应Y方向上的第round(l_Y,1·dec3))个位置。其中，round(·)表示取四舍五入的值，l_X,1表示第1帧观测数据的目标影响区域C₁中X方向的长度，l_Y,1表示第1帧观测数据的目标影响区域C₁中Y方向的长度。

(4.2)得出第1帧观测数据对应的每个初始粒子的权值，第1帧观测数据对应的第n个初始粒子的权值为：

其中，表示第1帧观测数据对应的第n个初始粒子值，表示第1帧观测数据的目标影响区域C₁中像素单元行数的集合，表示第1帧观测数据的目标影响区域C₁中像素单元列数的集合。为像素似然函数，其表达式为：

$l\left(z_1^{(i,j)}\right|{\hat{x}}_1^{\left(n\right)})=\exp [-\frac{h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)(h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)-2z_1^{(i,j)})}{{2\mathrm{\&sigma;}}_n^2}]$

其中，是观测噪声方差，表示坐标为(x₁,y₁)，强度为I₁的目标对位于第i行第j列的像素单元的强度贡献，其表达式为：

$h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)=\frac{\mathrm{\&Delta;x\&Delta;y}{I}_1}{2\mathrm{\&pi;}{\mathrm{\&Gamma;}}^2}\exp [-\frac{{(\mathrm{i\&Delta;x}-x_1)}^2+{(\mathrm{j\&Delta;y}-y_1)}^2}{{2\mathrm{\&Gamma;}}^2}]$

其中，Δx表示每帧观测数据对应的X方向分辨率，Δy表示每帧观测数据对应的Y方向分辨率。Γ为已知的雷达模糊斑点数量(参见“基于多模粒子滤波的机动弱目标检测前跟踪”中公式(10)中的Σ，龚亚信杨宏文胡卫东，2008,4)。

(4.3)对第1帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第1帧观测数据对应的每个初始粒子的归一化权值。第1帧观测数据对应的第n个初始粒子的归一化权值为：

$w_1^{\left(n\right)}=\frac{{\tilde{w}}_1^{\left(n\right)}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\tilde{w}}_1^{\left(n\right)}}$

(4.4)将第1帧观测数据对应的每个初始粒子按归一化权值大小进行排序，在第1帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第1帧观测数据对应的继续粒子，N_c为自然数且N_c<N。

(4.5)第1帧观测数据对应的每个继续粒子进行重采样，得出第1帧观测数据对应的N_c个重采样粒子、以及第1帧观测数据对应的每个重采样粒子的权值。第1帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第1帧观测数据对应的每个重采样粒子的权值相等，均为

然后，计算第1帧观测数据对应的目标状态值x₁：

$x_1=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_1^{\left(m\right)}{x}_1^{\left(m\right)}.$

步骤5，当k>1时，将第k-1帧观测数据对应的第m个重采样粒子代入重要性函数q(x_k|x_k-1,z_k)中进行采样，得到第k帧观测数据对应的第m个预测粒子其中，重要性函数q(x_k|x_k-1,z_k)选取先验分布p(x_k|x_k-1)。将m从1至N_c进行遍历，得出第k帧观测数据对应的N_c个预测粒子。按均匀分布密度函数抽取N_b个出生粒子，N_b＝N-N_c。这样便改善了粒子的多样性，充分利用了数据信息，保证目标的持续跟踪。

步骤6，得出第k帧观测数据对应的N个初始粒子，第k帧观测数据对应的N个初始粒子包括步骤5中第k帧观测数据对应的N_c个预测粒子和N_b个出生粒子；得出第k帧观测数据对应的每个初始粒子的权值；

其具体子步骤为：

得出第k帧观测数据对应的N个初始粒子，第k帧观测数据对应的N个初始粒子包括步骤5中第k帧观测数据对应的N_c个预测粒子和N_b个出生粒子。

得出第k帧观测数据对应的每个初始粒子的权值，第k帧观测数据对应的第n个初始粒子的权值为：

其中，表示第k帧观测数据对应的第n个初始粒子值，表示第k帧观测数据的目标影响区域C_k中像素单元行数的集合，表示第k帧观测数据的目标影响区域C_k中像素单元列数的集合。为像素似然函数，其表达式为：

$l\left(z_k^{(i,j)}\right|{\hat{x}}_k^{\left(n\right)})=\exp [-\frac{h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)(h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)-2z_k^{(i,j)})}{{2\mathrm{\&sigma;}}_n^2}]$

其中，是观测噪声方差，表示坐标为(x_k,y_k)，强度为I_k的目标对位于第i行第j列的像素单元的强度贡献，其表达式为：

$h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)=\frac{\mathrm{\&Delta;x\&Delta;y}{I}_k}{2\mathrm{\&pi;}{\mathrm{\&Gamma;}}^2}\exp [-\frac{{(\mathrm{i\&Delta;x}-x_k)}^2+{(\mathrm{j\&Delta;y}-y_k)}^2}{{2\mathrm{\&Gamma;}}^2}]$

其中，Δx表示每帧观测数据对应的X方向分辨率，Δy表示每帧观测数据对应的Y方向分辨率。Γ为已知的雷达模糊斑点数量(参见“基于多模粒子滤波的机动弱目标检测前跟踪”中公式(10)中的Σ，龚亚信杨宏文胡卫东，2008,4)。

步骤7，对第k帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第k帧观测数据对应的每个初始粒子的归一化权值。第k帧观测数据对应的第n个初始粒子的归一化权值为：

$w_k^{\left(n\right)}=\frac{{\tilde{w}}_k^{\left(n\right)}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\tilde{w}}_k^{\left(n\right)}}$

步骤8，将第k帧观测数据对应的每个初始粒子按归一化权值大小进行排序，在第k帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第k帧观测数据对应的继续粒子，N_c为自然数且N_c<N。

步骤9，对第k帧观测数据对应的每个继续粒子进行重采样，得出第k帧观测数据对应的N_c个重采样粒子、以及第k帧观测数据对应的每个重采样粒子的权值。第k帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第k帧观测数据对应的每个重采样粒子的权值相等，均为

然后，计算第k帧观测数据对应的目标状态值x_k：

$x_k=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_k^{\left(m\right)}{x}_k^{\left(m\right)}.$

步骤10，根据步骤5至步骤9，依次得出每帧观测数据对应的目标状态值，然后，根据每帧观测数据对应的目标状态值，输出目标航迹。

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

1)仿真实验环境

实验环境：MATLAB R2009b,Intel(R)Pentium(R)2CPU2.7GHz，Window7旗舰版。

2)仿真实验内容：利用某地面雷达回波仿真数据，应用本发明和传统标准粒子滤波方法进行目标检测跟踪。在仿真实验中，雷达脉冲重复频率为10KHz，脉冲宽度为10s，采样频率1MHz，带宽0.5MHz，扫描周期为1s。目标初始位置在距离雷达3km处，前10帧扫描时刻(对应前10帧观测数据)内目标以150m/s的速度在雷达径向30度方向上远离雷达匀速运动，后10帧扫描时刻(对应后10帧观测数据)内目标作变速曲线运动。在传统标准粒子滤波方法中，总粒子数为4000个，有效粒子数为500个。在本发明的改进粒子滤波方法中，初始粒子数为4000个(即本发明中的N)，继续粒子数为2000(即本发明中的N_c)个，出生粒子产生基数为2000(即本发明中的N_b)个，权值门限为欲淘汰粒子的平均值，恒虚警概率为10^-2。

将本发明和传统标准粒子滤波方法分别通过100次蒙特卡罗仿真平均，在不同信噪比环境下进行仿真。参照图2，为环境信噪比为3dB时两种方法得出的检测结果比较示意图，图2中，横轴表示水平距离，单位为米，纵轴表示垂直距离，单位为米。带方框的线条表示目标真实运动轨迹，带圆圈的线条表示应用传统标准粒子滤波方法得出的目标轨迹(对应图2图例中标准粒子滤波法检测轨迹)，带星号的线条表示应用本发明得出目标轨迹(对应图2图例中改进粒子滤波法检测轨迹)。参照图3，为环境信噪比为3dB时两种方法得出的距离检测均方误差比较示意图。图3中，横轴表示时间，单位为秒，纵轴表示距离，单位为米。图3中，带星号的线条表示应用传统标准粒子滤波方法得出的距离检测均方误差(对应图3图例中标准粒子滤波算法)，带方框线条表示应用本发明得出的距离检测均方误差(对应图3图例中改进粒子滤波算法)。

参照图4，为环境信噪比为5dB时两种方法得出的检测结果比较示意图，图4中，横轴表示水平距离，单位为米，纵轴表示垂直距离，单位为米。图4中，带方框的线条表示目标真实运动轨迹，带圆圈的线条表示应用传统标准粒子滤波方法得出的目标轨迹(对应图4图例中标准粒子滤波法检测轨迹)，带星号的线条表示应用本发明得出目标轨迹(对应图4图例中改进粒子滤波法检测轨迹)。参照图5，为环境信噪比为3dB时两种方法得出的距离检测均方误差比较示意图。图5中，横轴表示时间，单位为秒，纵轴表示距离，单位为米。图5中，带星号的线条表示应用传统标准粒子滤波方法得出的距离检测均方误差(对应图5图例中标准粒子滤波算法)，带方框线条表示应用本发明得出的距离检测均方误差(对应图5图例中改进粒子滤波算法)。

3)实验结果分析

在环境信噪比为3dB时，从图2可以看出，本发明相比于标准粒子滤波方法更加接近于真实的目标轨迹；从图3可以看出，在距离检测均方误差对比中可以看出，本发明相比于标准粒子滤波方法的距离检测均方误差更低。

在环境信噪比为5dB时，从图4可以看出，本发明相比于标准粒子滤波方法更加接近于真实的目标轨迹；从图5可以看出，在距离检测均方误差对比中可以看出，本发明相比于标准粒子滤波方法的距离检测均方误差更低。

综上所述，本发明首先通过基于自动删除算法的最小选择恒虚警率处理方法，对距离脉压后回波数据进行预处理，使粒子尽量分布在受目标影响的较小区域，减小搜索范围，从而减少了计算量；然后利用Halton拟随机序列在目标影响区域产生初始粒子，以代替标准粒子滤波方法中的蒙特卡洛随机采样，避免粒子产生群聚现象，使粒子分布的更加均匀，精度更高；最后将粒子划分为继续粒子和自适应出生粒子，保留权值较大的继续粒子进行动态更新，淘汰权值较小的粒子，同时根据粒子权值的大小分布情况，自适应地产生出生粒子，从而增加了粒子多样性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.基于改进粒子滤波的目标跟踪方法，其特征在于，包括以下步骤：

步骤1，当目标在一个平面内运动时，利用雷达接收目标的回波数据；对目标的回波数据进行距离向脉冲压缩处理，得到第1帧观测数据至第K帧观测数据，K为自然数；

步骤2，针对每帧观测数据，确定恒虚警检测门限；

步骤3，将每帧观测数据中低于恒虚警检测门限的数据丢弃，保留每帧观测数据中大于或等于恒虚警检测门限的数据；将第k帧观测数据中保留数据对应的区域记为第k帧观测数据的目标影响区域C_k，k取1至K；

步骤4，当k＝1时，在第k帧观测数据的目标影响区域C_k用Halton拟随机序列均匀地产生N个初始粒子，N为设定的自然数；得出第1帧观测数据对应的每个初始粒子的权值；对第1帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第1帧观测数据对应的每个初始粒子的归一化权值；在第1帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个初始粒子作为第1帧观测数据对应的继续粒子，N_c为自然数且N_c<N；对第1帧观测数据对应的每个继续粒子进行重采样，得出第1帧观测数据对应的N_c个重采样粒子、以及第1帧观测数据对应的每个重采样粒子的权值；第1帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第1帧观测数据对应的每个重采样粒子的权值为计算出第1帧观测数据对应的目标状态值x₁：

$x_1=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_1^{\left(m\right)}{x}_1^{\left(m\right)};$

步骤5，当k>1时，将第k-1帧观测数据对应的第m个重采样粒子代入重要性函数q(x_k|x_k-1,z_k)中进行采样，得到第k帧观测数据对应的第m个预测粒子将m从1至N_c进行遍历，得出第k帧观测数据对应的N_c个预测粒子；按均匀分布密度函数抽取N_b个出生粒子，N_b＝N-N_c；

步骤6，得出第k帧观测数据对应的N个初始粒子，第k帧观测数据对应的N个初始粒子包括步骤5中第k帧观测数据对应的N_c个预测粒子和N_b个出生粒子；得出第k帧观测数据对应的每个初始粒子的权值；

步骤7，对第k帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第k帧观测数据对应的每个初始粒子的归一化权值；

步骤8，在第k帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第k帧观测数据对应的继续粒子；

步骤9，对第k帧观测数据对应的每个继续粒子进行重采样，得出第k帧观测数据对应的N_c个重采样粒子、以及第k帧观测数据对应的每个重采样粒子的权值；第k帧观测数据对应的第m个重采样粒子表示为第k帧观测数据对应的每个重采样粒子的权值为计算出第k帧观测数据对应的目标状态值x_k：

$x_k=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_k^{\left(m\right)}{x}_k^{\left(m\right)}.$

2.如权利要求1所述的基于改进粒子滤波的目标跟踪方法，其特征在于，在步骤1中，以目标的运动平面为基准建立二维平面直角坐标系，二维平面直角坐标系的横轴为X轴，纵轴为Y轴；

第k帧观测数据表示为z_k， $z_k=\left\{z_k^{(i,j)}\right|1\&le;i\&le;M,1\&le;j\&le;N,1\&le;k\&le;K\},$ 表示第k帧观测数据中位于第i行第j列的像素单元对应的观测数据；M表示每帧观测数据对应的像素单元的行数，N表示每帧观测数据对应的像素单元的列数。

3.如权利要求2所述的基于改进粒子滤波的目标跟踪方法，其特征在于，所述步骤4的具体子步骤为：

(4.1)当k＝1时，在第k帧观测数据的目标影响区域C_k用Halton拟随机序列均匀地产生N个初始粒子的过程为：将十进制数dec1分别转化为2进制数和3进制数，dec1为自然数且dec1的取值为1至N；分别对2进制数和3进制数求反序，得到2进制反序后数和3进制反序后数；在2进制反序后数依次前置一个小数点和0，得到2进制小数，同时，在3进制反序后数依次前置一个小数点和0，得到3进制小数；将2进制小数转化为十进制数dec2，将3进制小数转化为十进制数dec3；则产生的第dec1个初始粒子对应X方向上的第round(l_X,1·dec2))个位置，同时对应Y方向上的第round(l_Y,1·dec3))个位置，round(·)表示取四舍五入，l_X,1表示第1帧观测数据的目标影响区域C₁中X方向的长度，l_Y,1表示第1帧观测数据的目标影响区域C₁中Y方向的长度；

(4.2)得出第1帧观测数据对应的每个初始粒子的权值，第1帧观测数据对应的第n个初始粒子的权值为：

其中，表示第1帧观测数据对应的第n个初始粒子值，表示第1帧观测数据的目标影响区域C₁中像素单元行数的集合，表示第1帧观测数据的目标影响区域C₁中像素单元列数的集合；为：

$l\left(z_1^{(i,j)}\right|{\hat{x}}_1^{\left(n\right)})=\exp [-\frac{h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)(h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)-2z_1^{(i,j)})}{{2\mathrm{\&sigma;}}_n^2}]$

其中，是观测噪声方差，的含义，表示坐标为(x₁,y₁)，强度为I₁的目标对位于第i行第j列的像素单元的强度贡献，其表达式为：

$h_1^{(i,j)}\left({\hat{x}}_1^{\left(n\right)}\right)=\frac{\mathrm{\&Delta;x\&Delta;y}{I}_1}{2\mathrm{\&pi;}{\mathrm{\&Gamma;}}^2}\exp [-\frac{{(\mathrm{i\&Delta;x}-x_1)}^2+{(\mathrm{j\&Delta;y}-y_1)}^2}{{2\mathrm{\&Gamma;}}^2}]$

其中，Δx表示每帧观测数据对应的X方向分辨率，Δy表示每帧观测数据对应的Y方向分辨率，Γ为已知的雷达模糊斑点数量；

(4.3)对第1帧观测数据对应的每个初始粒子的权值进行归一化处理，得到第1帧观测数据对应的每个初始粒子的归一化权值；第1帧观测数据对应的第n个初始粒子的归一化权值为：

$w_1^{\left(n\right)}=\frac{{\tilde{w}}_1^{\left(n\right)}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\tilde{w}}_1^{\left(n\right)}}$

其中，n取1至N；

(4.4)在第1帧观测数据对应的所有初始粒子中，抽取归一化权值最大的N_c个粒子作为第1帧观测数据对应的继续粒子；

(4.5)第1帧观测数据对应的每个继续粒子进行重采样，得出第1帧观测数据对应的N_c个重采样粒子、以及第1帧观测数据对应的每个重采样粒子的权值；第1帧观测数据对应的第m个重采样粒子表示为m＝1,…,N_c，第1帧观测数据对应的每个重采样粒子的权值为计算出第1帧观测数据对应的目标状态值x₁：

$x_1=\underset{m=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{w}_1^{\left(m\right)}{x}_1^{\left(m\right)}.$

4.如权利要求2所述的基于改进粒子滤波的目标跟踪方法，其特征在于，在步骤6中，第k帧观测数据对应的第n个初始粒子的权值为：

其中，表示第k帧观测数据对应的第n个初始粒子值，表示第k帧观测数据的目标影响区域C_k中像素单元行数的集合，表示第k帧观测数据的目标影响区域C_k中像素单元列数的集合；为为：

$l\left(z_k^{(i,j)}\right|{\hat{x}}_k^{\left(n\right)})=\exp [-\frac{h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)(h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)-2z_k^{(i,j)})}{{2\mathrm{\&sigma;}}_n^2}]$

其中，是观测噪声方差，表示坐标为强度为I_k的目标对位于第i行第j列的像素单元的强度贡献，其表达式为：

$h_k^{(i,j)}\left({\hat{x}}_k^{\left(n\right)}\right)=\frac{\mathrm{\&Delta;x\&Delta;y}{I}_k}{2\mathrm{\&pi;}{\mathrm{\&Gamma;}}^2}\exp [-\frac{{(\mathrm{i\&Delta;x}-x_k)}^2+{(\mathrm{j\&Delta;y}-y_k)}^2}{{2\mathrm{\&Gamma;}}^2}]$

其中，Δx表示每帧观测数据对应的X方向分辨率，Δy表示每帧观测数据对应的Y方向分辨率；Γ为已知的雷达模糊斑点数量。

5.如权利要求1所述的基于改进粒子滤波的目标跟踪方法，其特征在于，在步骤7中，第k帧观测数据对应的第n个初始粒子的归一化权值为：

$w_k^{\left(n\right)}=\frac{{\tilde{w}}_k^{\left(n\right)}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\tilde{w}}_k^{\left(n\right)}}$

其中，n取1至N，为第k帧观测数据对应的第n个初始粒子的权值。

6.如权利要求1所述的基于改进粒子滤波的目标跟踪方法，其特征在于，在步骤9之后，根据每帧观测数据对应的目标状态值，得出目标航迹。