CN104237853A - 一种针对多帧检测前跟踪点迹序列的粒子滤波方法 - Google Patents

Abstract

本发明提供一种针对多帧检测前跟踪点迹序列的粒子滤波方法。本法米国内发明利用多帧检测前跟踪处理得到的短航迹对目标状态进行滤波处理，可以提高了对目标状态估计的精度；对多帧检测前跟踪处理得到的短航迹做进一步处理，实现了对目标的长时间跟踪，解决了多帧检测前跟踪方法不能提供目标完整航迹信息的问题；可以迭代的对目标状态进行实时估计，并且没有引入过多的存储量和计算量；通过粒子状态预测，可以有效解决低信噪比微弱目标检测前跟踪可能出现的漏检问题，保证了航迹的连续性。

Description

一种针对多帧检测前跟踪点迹序列的粒子滤波方法

技术领域

本发明属于雷达数据处理技术领域，它特别涉及了雷达微弱目标跟踪技术领域。

背景技术

雷达探测环境的日益复杂和目标隐身技术的日益成熟，使得雷达及时可靠地监测到淹没在强杂波中的微弱目标变得更加困难。提高雷达对微弱目标的检测跟踪性能具有重要的理论价值和实际意义。

检测前跟踪技术是一种在低信噪比环境中对微弱目标进行检测跟踪的一种技术。与传统的检测方法不同，检测前跟踪不单帧宣布量测结果，而是对多帧数据进行处理后同时宣布目标的量测结果与航迹。实现检测前跟踪技术的方法包括有粒子滤波、动态规划、霍夫变换等。基于粒子滤波的检测前跟踪方法的航迹起始能力不如动态规划和霍夫变换等方法，并且该方法需要计算概率，计算量大。但是，动态规划等多帧检测前跟踪方法处理一次只能得到目标的短航迹，无法得到目标的完整航迹。而且，一些多帧检测前跟踪方法是在离散空间得到的状态估计结果，估计精度存在一定损失。

发明内容

本发明所要解决的技术问题是，为提高对目标状态的估计精度，提供一种用于处理多帧检测前跟踪得到的短航迹的滤波方法。

本发明为解决上述技术问题所采用的技术方案是，一种针对多帧检测前跟踪点迹序列的粒子滤波方法，包括以下步骤：

步骤1)变量初始化：

初始化总仿真时间K、一次检测前跟踪处理帧数N；检测前跟踪处理门限V_T、检测前跟踪量测可能范围γ、一次多帧检测前跟踪处理第n帧的检测概率n＝1,2,…,N、雷达扫描周期T、目标状态转移矩阵F、观测矩阵H、过程噪声协方差Q、粒子数N_s、x和y方向的距离范围[x_min,x_max]和[y_min,y_max]、x和y方向的速度范围和

步骤2)粒子集初始化：

初始化时间变量k＝1与第i个粒子状态在x和y方向位置和速度状态 $x_k^i={[x_k^i,{\stackrel{\·}{x}}_k^i,y_k^i,{\stackrel{\·}{y}}_k^i]}^T,i=\mathrm{1,2},...,N_s{[\·]}^T$ 表示转置；

初始化粒子集初始化完成后，进入步骤4)；

步骤3)状态预测：

利用k-1时刻的粒子集对k时刻进行预测：得到预测粒子集 是过程噪声协方差Q的过程噪声，F为目标状态转移矩阵；

当第i个粒子状态中元素大于其最大范围的，则将最大范围值赋值给该元素，并设置该粒子对应的权值置为0；

步骤4)多帧检测前跟踪处理：

从雷达接收机中读取第k-N+1,k-N+2,…,k帧回波数据，当k≤N，则从雷达接收机中读取第1,2,…,k帧回波数据；当第k次检测前跟踪处理中存在量测结果，则进入步骤5)，否则进入步骤7)；存在量测结果即表示量测结果超过门限；

步骤5)权值更新：

k≤N，则量测结果为 $z_{k-N+1,k-N+2,...,k}^{t_k}=[z_{k-N+1}^{t_k},z_{k-N+2}^{t_k},...,z_k^{t_k}],$ 则更新粒子权值：

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s;$

当k>N，量测结果为开始更新粒子权值：

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_k-N+1})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_k-N+2})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s;$

表示最后一帧时刻为t_k的检测前跟踪处理得到的第l帧量测结果，表示对应父粒子的状态，当l≤k-1，用表示对应父粒子的状态；P表示

似然概率；

表示向下取整，H为观测矩阵，量测范围变量△_x和△_y的取值范围均为且△_x和△_y均不等于零；当γ＝1时似然概率的计算式中无中间项；

计算更新后的粒子权值和再对粒子权值归一化处理进入步骤6)

步骤6)对权值更新后的粒子集进行***重采样；

6.1)初始化

6.2)对于所有i＝2,3,…,N_s，计算

6.3)产生随机数u₁，u₁服从范围内的均匀分布；

6.4)令i＝1，对于所有j＝1,2,…,N_s，执行计算：

当u_j>c_i，则更新i＝i+1，返回步骤6.4)直至u_j≤c_i；再令并记录该粒子的父粒子i^j＝i；

重采样后的粒子集表示为进入步骤7)

步骤7)目标状态估计：

对所有重采样后的粒子状态加权平均得到第k时刻的目标状态估计结果；

步骤8)更新k＝k+1，当k≤K，返回步骤3)，否则N帧粒子滤波结束。

本发明通过上面的步骤，就可以根据检测前跟踪处理得到的短航迹更精确地实时估计出目标状态，得到目标的完整航迹。

本发明的有益效果是，本发明利用多帧检测前跟踪处理得到的短航迹对目标状态进行滤波处理，可以提高了对目标状态估计的精度；对多帧检测前跟踪处理得到的短航迹做进一步处理，实现了对目标的长时间跟踪，解决了多帧检测前跟踪方法不能提供目标完整航迹信息的问题；可以迭代的对目标状态进行实时估计，并且没有引入过多的存储量和计算量；通过粒子状态预测，可以有效解决低信噪比微弱目标检测前跟踪可能出现的漏检问题，保证了航迹的连续性。

附图说明

图1为本发明的流程图。

图2为本发明实施例中目标真实航迹与应用本发明跟踪得到的目标航迹。

图3为本发明实施例中应用本发明方法估计目标状态的目标位置均方根误差。

图4为本发明实施例中应用本发明方法估计目标状态的目标速度均方根误差。

具体实施方式

如图1所示，本发明流程如下：

步骤1)变量初始化：

初始化总仿真时间K；一次检测前跟踪处理帧数N；目标信噪比SNR；检测前跟踪处理门限VT；检测前跟踪量测可能范围γ，γ∈{δ²；δ＝1,2,3,4,…}，δ表示检测单元数，检测前跟踪量测可能范围表示为检测前跟踪的量测结果与真实目标之间不能超过的范围；一次多帧检测前跟踪处理第n帧的检测概率n＝1,2,…,N，检测概率在进行检测前跟踪时已知；雷达扫描周期T；目标状态转移矩阵F；观测矩阵H；过程噪声协方差Q；粒子数N_s；x和y方向的距离范围[x_min,x_max]和[y_min,y_max]；x和y方向的速度范围和

步骤2)粒子集初始化：

初始化时间变量k＝1与第i个粒子状态在x和y方向位置和速度状态

$x_k^i=x_{\min }+(x_{\max }-x_{\min })\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

${\stackrel{\·}{x}}_k^i={\stackrel{\·}{x}}_{\min }+({\stackrel{\·}{x}}_{\max }-{\stackrel{\·}{x}}_{\min })\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

$y_k^i=y_{\min }+(y_{\max }-y_{\min })\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

${\stackrel{\·}{y}}_k^i={\stackrel{\·}{y}}_{\min }+({\stackrel{\·}{y}}_{\max }-{\stackrel{\·}{y}}_{\min })\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

rand(1)表示按照均匀分布产生一个[0,1]区间上的随机数；

第i个粒子权重初始化为： $w_k^i=\frac{1}{N_s},i=\mathrm{1,2},...,N_s;$

粒子集初始化完成后，进入步骤4)；

步骤3)状态预测：

利用k-1时刻的粒子集对k时刻进行预测：得到预测粒子集 是过程噪声协方差Q的过程噪声；

当第i个粒子状态中元素大于其最大范围的，则将最大范围值赋值给该元素，并设置该粒子对应的权值置为0；

步骤4)多帧检测前跟踪处理：

从雷达接收机中读取第k-N+1,k-N+2,…,k帧回波数据，当k≤N，则从雷达接收机中读取第1,2,…,k帧回波数据；当第k次检测前跟踪处理中存在量测结果，则进入步骤5)，否则进入步骤7)；

步骤5)权值更新：

k≤N，则量测结果为 $z_{k-N+1,k-N+2,...,k}^{t_k}=[z_{k-N+1}^{t_k},z_{k-N+2}^{t_k},...,z_k^{t_k}],$ 则更新粒子权值：

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s$

当k>N，量测结果为开始更新粒子权值：

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_k-N+1})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_k-N+2})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s$

表示最后一帧时刻为t_k的检测前跟踪处理得到的第l帧量测结果，表示对应父粒子的状态，当l≤k-1，用表示对应父粒子的状态；P表示

似然概率；

表示向下取整，H为观测矩阵，量测范围变量△_x和△_y的取值范围均为且△_x和△_y均不等于零；

计算更新后的粒子权值和再对粒子权值归一化处理进入步骤6)

步骤6)对权值更新后的粒子集进行***重采样；

6.1)初始化

6.2)对于所有i＝2,3,…,N_s，计算

6.3)产生随机数u₁，u₁服从范围内的均匀分布；

6.4)令i＝1，对于所有j＝1,2,…,N_s，执行计算：

当u_j>c_i，则更新i＝i+1，返回步骤6.4)直至u_j≤c_i；再令并记录该粒子的父粒子i^j＝i；

重采样后的粒子集表示为进入步骤7)

步骤7)目标状态估计：

对所有重采样后的粒子状态加权平均得到第k时刻的目标状态估计结果；

步骤8、更新k＝k+1，当k≤K，返回步骤3)，否则N帧粒子滤波结束。

实施例

本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2010b上验证正确。具体实施步骤如下：

步骤1、变量初始化：

初始化变量包括：总仿真时间K＝30；一次检测前跟踪处理帧数N＝6；目标信噪比SNR＝8dB；检测前跟踪处理门限V_T＝13.2847；检测前跟踪量测可能范围γ＝7²；一次多帧检测前跟踪处理第第1,2,…,N帧的检测概率分别为0.876,0.892,0.918,0.910,0.908,0.838；雷达扫描周期T＝1s；目标状态转移矩阵F

$F=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]$

观测矩阵H

$H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

过程噪声协方差Q

$Q=\left[\begin{array}{cccc}0.001\×T^3/3& 0.001\×T^2/2& 0& 0\\ 0.001\×T^2/2& 0.001\×T& 0& 0\\ 0& 0& 0.001\×T^3/3& 0.001\×T^2/2\\ 0& 0& 0.001\×T^2/2& 0.001\×T\end{array}\right]$

粒子数N_s＝8000；x和y方向的距离范围均为[1,59]；x和y方向的速度范围均为[0.6,1.2]。

步骤2、粒子集初始化：

初始化时间变量k＝1。粒子x和y方向位置和速度状态初始化

$x_k^i=1+58\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

${\stackrel{\·}{x}}_k^i=0.6+0.6\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

$y_k^i=1+58\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

${\stackrel{\·}{y}}_k^i=0.6+0.6\×\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

其中rand(1)表示按照均匀分布产生一个[0,1]区间上的随机数。令

$x_k^i={[x_k^i,{\stackrel{\·}{x}}_k^i,y_k^i,{\stackrel{\·}{y}}_k^i]}^T,i=\mathrm{1,2},...,N_s$

粒子权重初始化

$w_k^i=\frac{1}{N_s},i=\mathrm{1,2},...,N_s$

则，得到初始粒子集

步骤3、状态预测：

如果k＝1，则跳过此步骤。否则对k-1时刻的粒子集进行一步预测

$x_{k|k-1}^i={\mathrm{Fx}}_{k-1}^i+v_k^i,i=\mathrm{1,2},...,N_s$

得到预测粒子集其中是均值为零协方差矩阵为Q的过程噪声。

如果第i个粒子状态 $x_{k|k-1}^i\left(1\right)<1$ 或者 $x_{k|k-1}^i\left(1\right)>59,$ 则令 $x_{k|k-1}^i\left(1\right)=1$ 或者 $x_{k|k-1}^i\left(1\right)=59,$ 如果 $x_{k|k-1}^i\left(3\right)<1$ 或者 $x_{k|k-1}^i\left(3\right)>59,$ 则令 $x_{k|k-1}^i\left(3\right)=1$ 或者 $x_{k|k-1}^i\left(3\right)=59.$ 并将该粒子对应的权值置为0。

步骤4、多帧检测前跟踪处理：

从雷达接收机中读取第k-N+1,k-N+2,…,k帧回波数据(如果k≤N，则从雷达接收机中读取第1,2,…,N帧回波数据)，进行多帧检测前跟踪处理。

步骤5、权值更新：

令表示最后一帧时刻为t_k的检测前跟踪处理得到的第l帧量测结果。如果在该次检测前跟踪处理中存在量测结果(如果k≤N，则量测结果为)，则更新预测粒子权值：

为了表示方便，用表示对应父粒子的状态，当l≤k-1，用表示对应父粒子的状态。当k>N时，按下式更新粒子权值

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_k-N+1})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_k-N+2})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s$

其中，似然概率

△_x和△_y的取值范围均为且△_x和△_y均不等于零。

当k≤N时，按下式更新粒子权值

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s$

计算粒子权值和

$\mathrm{\&Xi;}=\underset{i=1}{\overset{N_s}{\mathrm{\&Sigma;}}}{w}_k^i$

在Ξ不为零的情况下对粒子权值归一化

$w_k^i=\frac{w_k^i}{\mathrm{\&Xi;}},i=\mathrm{1,2},...,N_s$

如果在该次检测前跟踪处理中不存在量测结果，即没有状态过门限，则跳过该步骤。

步骤6、重采样：

如果在该次检测前跟踪处理中不存在量测结果，即没有状态过门限，则跳过该步骤。否则对权值更新后的粒子集进行***重采样。

6.1、初始化

6.2、对于所有i＝2,3,…,N_s，计算

$c_i=c_{i-1}+w_k^i$

6.3、随机产生u₁，其服从范围内的均匀分布，即

6.4、令i＝1，对于所有j＝1,2,…,N_s，执行以下计算：

$u_j=u_1+(j-1)N_s^{-1}$

如果u_j>c_i，则循环执行i:＝i+1，直至u_j≤c_i。然后令并记录该粒子的父粒子i^j＝i。

重采样后的粒子集表示为

步骤7、目标状态估计：

对所有粒子状态加权平均得到第k时刻的目标状态估计结果。

步骤8、令k:＝k+1，如果k≤K，重复步骤3～步骤7。

附图2给出了本仿真实例中对目标航迹的恢复结果，成功提供了目标的完整航迹信息，且不存在漏检问题。附图3给出了本实施例中利用本发明得到的目标位置均方根误差，附图4给出了本实施例中利用本发明得到的目标速度均方根误差，是200次蒙特卡罗试验的统计结果，结果表明该方法具有较高且平稳的估计精度。

Claims (2)

1.一种针对多帧检测前跟踪点迹序列的粒子滤波方法，其特征在于，包括以下步骤：

步骤1)变量初始化：

初始化总仿真时间K、一次检测前跟踪处理帧数N；检测前跟踪处理门限V_T、检测前跟踪量测可能范围γ、一次多帧检测前跟踪处理第n帧的检测概率n＝1,2,…,N、雷达扫描周期T、目标状态转移矩阵F、观测矩阵H、过程噪声协方差Q、粒子数N_s、x和y方向的距离范围[x_min,x_max]和[y_min,y_max]、x和y方向的速度范围和

步骤2)粒子集初始化：

初始化时间变量k＝1与第i个粒子状态在x和y方向位置和速度状态 $x_k^i={[x_k^i,{\stackrel{\&CenterDot;}{x}}_k^i,y_k^i,{\stackrel{\&CenterDot;}{y}}_k^i]}^T,i=\mathrm{1,2},...,N_s{[\&CenterDot;]}^T$ 表示转置；

初始化粒子集初始化完成后，进入步骤4)；

步骤3)状态预测：

利用k-1时刻的粒子集对k时刻进行预测：得到预测粒子集 是过程噪声协方差Q的过程噪声，F为目标状态转移矩阵；

当第i个粒子状态中元素大于其最大范围的，则将最大范围值赋值给该元素，并设置该粒子对应的权值置为0；

步骤4)多帧检测前跟踪处理：

从雷达接收机中读取第k-N+1,k-N+2,…,k帧回波数据，当k≤N，则从雷达接收机中读取第1,2,…,k帧回波数据；当第k次检测前跟踪处理中存在量测结果，则进入步骤5)，否则进入步骤7)；存在量测结果即表示量测结果超过门限；

步骤5)权值更新：

k≤N，则量测结果为 $z_{k-N+1,k-N+2,...,k}^{t_k}=[z_{k-N+1}^{t_k},z_{k-N+2}^{t_k},...,z_k^{t_k}],$ 则更新粒子权值：

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s;$

当k>N，量测结果为开始更新粒子权值：

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_k-N+1})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_k-N+2})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=\mathrm{1,2},...,N_s;$

表示最后一帧时刻为t_k的检测前跟踪处理得到的第l帧量测结果，表示对应父粒子的状态，当l≤k-1，用表示对应父粒子的状态；P表示

似然概率；

表示向下取整，H为观测矩阵，量测范围变量△_x和△_y的取值范围均为且△_x和△_y均不等于零；当γ＝1时似然概率的计算式中无中间项；

计算更新后的粒子权值和再对粒子权值归一化处理进入步骤6)

步骤6)对权值更新后的粒子集进行***重采样；

6.1)初始化

6.2)对于所有i＝2,3,…,N_s，计算

6.3)产生随机数u₁，u₁服从范围内的均匀分布；

6.4)令i＝1，对于所有j＝1,2,…,N_s，执行计算：

当u_j>c_i，则更新i＝i+1，返回判断u_j>c_i是否成立，如果成立则继续更新i＝i+1，直至u_j≤c_i；再令并记录该粒子的父粒子i^j＝i；

重采样后的粒子集表示为进入步骤7)

步骤7)目标状态估计：

对所有重采样后的粒子状态加权平均得到第k时刻的目标状态估计结果；

步骤8)更新k＝k+1，当k≤K，返回步骤3)，否则N帧粒子滤波结束。

2.如权利要求1所述一种针对多帧检测前跟踪点迹序列的粒子滤波方法，其特征在于，初始化时间变量k＝1与第i个粒子状态在x和y方向位置和速度状态的具体方法为：

$x_k^i=x_{\min }+(x_{\max }-x_{\min })\&times;\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

${\stackrel{\&CenterDot;}{x}}_k^i={\stackrel{\&CenterDot;}{x}}_{\min }+({\stackrel{\&CenterDot;}{x}}_{\max }-{\stackrel{\&CenterDot;}{x}}_{\min })\&times;\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

$y_k^i=y_{\min }+(y_{\max }-y_{\min })\&times;\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

${\stackrel{\&CenterDot;}{y}}_k^i={\stackrel{\&CenterDot;}{y}}_{\min }+({\stackrel{\&CenterDot;}{y}}_{\max }-{\stackrel{\&CenterDot;}{y}}_{\min })\&times;\mathrm{rand}\left(1\right),i=\mathrm{1,2},...,N_s$

rand(1)表示按照均匀分布产生一个[0,1]区间上的随机数；

初始化第i个粒子权重的具体方法为为：