CN102621542A - 基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法 - Google Patents

Abstract

本发明公开了一种基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法，属于雷达数据处理领域。本发明提出的方法立足于解决低信噪比情况下的机动微弱目标检测与跟踪问题，在实现对机动微弱目标检测同时提供目标航迹信息，并有效消除目标信噪比低、机动大造成的目标漏检问题。本发明的方法主要包括以下步骤：(1)利用多模粒子滤波得到各时刻的目标状态和状态协方差估计；(2)将多模粒子滤波的结果重新定义为量测值，利用“状态估计-航迹”数据关联给出目标航迹信息。本发明克服了基于多模粒子滤波的机动微弱目标检测前跟踪方法的局限性，通过有效降低漏检概率，进而保证目标航迹的连续性，同时本发明具有结构简单、易于硬件实现等优点，因此具有较强的工程应用价值和推广前景。

Description

基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法

技术领域

本发明涉及一种雷达数据处理方法，特别涉及一种雷达对机动微弱目标的检测与跟踪方法。

背景技术

微弱目标的检测与跟踪是现代雷达面临的难点问题之一。各种隐身飞机、掠海导弹以及巡航导弹等微弱目标的大量出现和使用，严重的影响着雷达目标检测跟踪性能的发挥。同时，为了提高生存和突防能力，微弱目标会随时进行一定的机动，其运动状态也会因不同的作战任务而变，进一步增大了雷达对其进行检测与跟踪的困难。因此，实现雷达对机动微弱目标的有效检测与跟踪，对提高雷达作战效能和生存能力具有重要意义。目前，基于粒子滤波的微弱目标检测前跟踪方法由于具有处理方法简单、易于硬件实现、可适应于非线性非高斯系统等优点，在微弱目标检测与跟踪领域受到广泛重视，而基于粒子滤波的检测前跟踪方法通过借鉴交互式多模型的思想，形成了一种基于多模粒子滤波的机动微弱目标检测前跟踪方法。该方法主要通过以下步骤实现：

(1)将雷达接收机输出的目标回波信号进行A/D变换，送雷达数据处理计算机执行以下步骤；

(2)对目标运动模式进行预测；

(3)利用目标运动模式和相应的目标状态转移方程对粒子集进行预测；

(4)利用新量测对各粒子权重进行更新；

(5)对权重更新后的粒子集进行重采样，并得到目标状态估计。

基于多模粒子滤波的机动微弱目标检测前跟踪方法具有以下两个缺陷：(1)不能提供目标的航迹信息；(2)当目标信噪比较低或机动较大时容易出现目标漏检。

发明内容

本发明的目的是提出一种基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法，解决基于多模粒子滤波的机动微弱目标检测前跟踪方法容易出现漏检和无法提供目标航迹信息的问题。

本发明提出的基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法的技术方案包括以下步骤：

步骤1：变量初始化

K是总仿真时间；

T是雷达扫描周期；

N是滤波器采用的粒子数；

SNR是目标信噪比；

I_min和I_max是目标可能的最小和最大强度；

v_min和v_max是目标可能的最小和最大速度；

γ是目标存在判决门限；

是目标初始运动模式

η₀是目标初始存在概率；

L×H表示雷达监测区域的大小；

∏_e是目标存在状态概率转移矩阵；

∏_m是目标运动模式概率转移矩阵；

Q_k是过程噪声协方差；

步骤2：令k＝0，进行粒子初始化，得到初始粒子集

其中，

包含目标的位置

速度和强度

信息，

是目标运动模式变量，是目标存在状态变量，i＝1，2，…，N

(1)令i＝1；

(2)目标位置初始化

$x_k^l=L\×\mathrm{rand}\left(1\right)---\left(1\right)$

$y_k^l=H\×\mathrm{rand}\left(1\right)---\left(2\right)$

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

(3)目标速度初始化

${\stackrel{\·}{x}}_k^l=v_{\min }+(v_{\max }-v_{\min })\×\mathrm{rand}\left(1\right)---\left(3\right)$

${\stackrel{\·}{y}}_k^l=v_{\min }+(v_{\max }-v_{\min })\×\mathrm{rand}\left(1\right)---\left(4\right)$

(4)目标强度初始化

$I_k^l=I_{\min }+(I_{\max }-I_{\min })\×\mathrm{rand}\left(1\right)---\left(5\right)$

(5)根据目标初始运动模式概率初始化目标运动模式变量

(6)根据目标初始存在概率η₀初始化目标存在状态变量

(7)粒子权重初始化

$W_k^l=\frac{1}{N}---\left(6\right)$

(8)令i＝i+1，若i≤N，转(2)；

步骤3：令k＝k+1，获得k时刻的雷达量测

将雷达接收到的信号进行A/D变换，得到k时刻的雷达量 $z_k=\left\{z_k^{(p,q)}\right|p=\mathrm{1,2},...,m;q=\mathrm{1,2},...,n\},$ 送雷达数据处理计算机；

步骤4：对初始粒子集进行预测，得到粒子集

其中，

(1)令i＝1；

(2)根据前一时刻的目标存在状态变量

和目标存在状态概率转移矩阵∏_e预测当前时刻目标存在状态变量

(3)若且

利用方程(1)-方程(5)生成粒子状态

并根据

随机生成目标运动模式变量

(4)若

且根据前一时刻的目标运动模式变量

和目标运动模式概率转移矩阵∏_m预测当前时刻目标运动模式变量

(5)根据

选择目标状态状态转移矩阵，其中

$F_k^{\left(1\right)}=\left[\begin{array}{ccccc}1& T& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& T& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]---\left(7\right)$

对应于目标匀速直线运动(CV模型)

$F_k^{\left(j\right)}=\left[\begin{array}{ccccc}1& \frac{\sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 0& -\frac{1-\cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 0\\ 0& \cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0& -\sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0\\ 0& \frac{1-\cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 1& \frac{\sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 0\\ 0& \sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0& \cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0\\ 0& 0& 0& 0& 1\end{array}\right],$ j＝2，3(8)

对应于目标坐标转动(CT模型)，其中

${\mathrm{\Ω}}_k^{\left(2\right)}=\frac{{-\mathrm{\α}}_m}{\sqrt{{\left({\stackrel{\·}{x}}_{k|k-1}^l\right)}^2+{\left({\stackrel{\·}{y}}_{k|k-1}^l\right)}^2}},$ ${\mathrm{\Ω}}_k^{\left(3\right)}=\frac{{\mathrm{\α}}_m}{\sqrt{{\left({\stackrel{\·}{x}}_{k|k-1}^l\right)}^2+{\left({\stackrel{\·}{y}}_{k|k-1}^l\right)}^2}}---\left(9\right)$

分别对应于顺时针坐标转动和逆时针坐标转动，α_m＞0是一个典型的机动加速度；

(6)根据目标运动模式变量

选择相应的状态转移方程，对粒子i的状态进行一步预测

$x_{k|l-1}^l=F_k^{\left(r_{k|k-1}\right)}{x}_{k-1}^l+v_k---\left(10\right)$

其中，v_k为过程噪声，其噪声协方差为Q_k；

(7)令i＝i+1，若i≤N，转(2)；

步骤5：利用当前时刻的量测对粒子权重进行更新

(1)令i＝1；

(2)确定似然区域的范围

$C_p\left(x_{k|k-1}^l\right)=\{p_0-s,...,p_0-1,p_0,p_0+1,...,p_0+s\}---\left(11\right)$

$C_q\left(y_{k|k-1}^l\right)=\{q_0-s,...,q_0,q_0+1,...,q_0+s\}---\left(12\right)$

其中，p₀表示最接近状态矢量元素的整数，q₀表示最接近状态矢量元素

的整数，s为预设参数；

(3)粒子权重更新

$w_{k|k-1}^l=\left\{\begin{array}{cc}\underset{p\∈C_p\left(x_{k|k-1}^l\right)}{\mathrm{\Π}}\underset{q\∈C_q\left(x_{k|k-1}^l\right)}{\mathrm{\Π}}l\left(z_k^{(p,q)}\right|x_{k|k-1}^l)& E_{k|k-1}^l=1\\ 1& E_{k|k-1}^l=0\end{array}\right.---\left(13\right)$

其中

$l\left(z_k^{(p,q)}\right|x_{k|k-1}^l)=\exp \{-\frac{h_k^{(p,q)}\left(x_{k|k-1}^l\right)(h_k^{(p,q)}\left(x_{k|k-1}^l\right)-{2z}_k^{(p,q)})}{{2\mathrm{\σ}}^2}\}---\left(14\right)$

$h_k^{(p,q)}\left(x_{k|k-1}^l\right)=\frac{{\mathrm{\Δ}}_x{\mathrm{\Δ}}_y{I}_{k|k-1}^l}{2\mathrm{\π}{\mathrm{\Σ}}^2}\exp \{-\frac{{({\mathrm{p\Δ}}_x-x_{k|k-1}^l)}^2+{({\mathrm{q\Δ}}_y-y_{k|k-1}^l)}^2}{2{\mathrm{\Σ}}^2}\}---\left(15\right)$

其中，σ²为量测噪声协方差，Δ_x和Δ_y表示传感器分辨率，∑表示传感器模糊斑点数；

(4)令i＝i+1，若i≤N，转(2)；

(5)权重归一化

$w_k^l=\frac{w_{k|k-1}^l}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{k|k-1}^j},$ i＝1，2，…，N    (16)

步骤6：重采样

对权重更新后的粒子集进行重采样

${\{y_k^l,(1/N)\}}_{l=1}^N=\mathrm{resample}{\left\{y_{k|k-1}^l{w}_k^l\right\}}_{l=1}^N---\left(17\right)$

得到k时刻的粒子集其中

步骤7：目标状态和状态协方差估计

(1)计算目标存在概率

${\hat{P}}_k=\frac{1}{N}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{E}_k^l---\left(16\right)$

(2)若进行目标状态估计

${\hat{x}}_k=\frac{\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^l{E}_k^l}{\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{E}_k^l}---\left(19\right)$

和状态协方差估计

${\hat{B}}_k=\frac{1}{N}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}({\hat{x}}_k-x_k^l){({\hat{x}}_k-x_k^l)}^T---\left(20\right)$

步骤8：“状态估计-航迹”数据关联

(1)如果当前时刻不存在目标状态估计

转(4)；

(2)将目标状态估计

等效为新的量测值，状态协方差估计

等效为量测误差；

(3)若当前时刻不存在目标航迹，将新量测作为航迹头，利用m/n逻辑法进行航迹起始；否则，利用新量测对目标航迹τ_k-1或τ′_k-1进行更新，得到更新的航迹τ_k，转步骤9；

(4)如果当前时刻存在目标航迹τ_k-1或τ′_k-1，若已经连续3个时刻不存在目标状态估计，输出航迹τ_k-3并删除航迹τ′_k-1；否则，对τ_k-1或τ′_k-1进行一步预测，得到预测τ′_k；

步骤9：重复步骤3～步骤8，直至雷达关机。

和背景技术相比，本发明的有益效果说明：

(1)本发明提出的基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法，可以在对目标状态估计的同时给出目标航迹，解决了背景技术中基于多模粒子滤波的机动微弱目标检测前跟踪方法不能保持目标航迹信息的问题；(2)本发明提出的方法通过对“状态估计-航迹”数据关联，可以有效消除目标信噪比低、机动大造成的某些时刻的目标漏检问题，进而保证目标航迹的连续性。

附图说明

附图1是本发明的基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法的整体流程图，附图中各符号的含义与发明内容部分相应符号的含义相同；

附图2是本发明实施例中基于多模粒子滤波的机动微弱目标检测前跟踪方法估计的目标存在概率；

附图3是本发明实施例中目标真实航迹与基于多模粒子滤波的机动微弱目标检测前跟踪方法估计的目标状态，附图中X位置和Y位置分别表示目标的x坐标和y坐标；

附图4是本发明实施例中目标真实航迹与基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法估计的目标航迹，附图中X位置和Y位置分别表示目标的x坐标和y坐标；

具体实施方式

下面结合附图对本发明的基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法进行详细描述。

不失一般性，设置一个二维仿真场景。假设目标初始存在概率η₀＝0.05，目标初始强度为I＝20，其可能的最小强度和最大强度分别为I_min＝10和I_max＝30，目标初始位置和初始速度分别为(20，20)和(0.45，0.25)，可能的最小速度和最大速度分别为v_min＝-1和v_max＝1，目标状态转移概率P_b＝P_d＝0.05，目标存在判决门限γ＝0.6，目标从第7s开始出现，首先沿x轴方向作10s匀速直线运动，接着作5s逆时针转弯运动，此后进入10s的匀速直线运动，最后作5s的顺时针机动转弯，直到37s消失，过程噪声强度分别为q₁＝0.001和q₂＝0.01；雷达扫描周期T＝1，分辨率Δx＝Δy＝1，分辨单元的格式m＝n＝60，量测噪声标准差σ＝2.90，传感器模糊斑点数∑＝0.7；滤波器采用的粒子数为N＝40000，总仿真时间为K＝45个扫描周期，目标似然区域选择参数s＝2。本发明具体步骤如附图1所示。

(1)根据以上仿真条件进行变量初始化：

步骤1：变量初始化

K＝45；

T＝1；

N＝4000；

SNR＝7dB；

I_min＝10，I_max＝30；

v_min＝-2，v_max＝2；

γ＝0.6；

η₀＝0.05；

L＝60，H＝60；

${\mathrm{\Π}}_e=\left[\begin{array}{cc}0.95& 0.05\\ 0.05& 0.95\end{array}\right];$

${\mathrm{\Π}}_m=\left[\begin{array}{ccc}0.9& 0.05& 0.05\\ 0.4& 0.5& 0.1\\ 0.4& 0.1& 0.5\end{array}\right];$

$Q_k=\left[\begin{array}{ccc}{Q}_0& 0& 0\\ 0& Q_0& 0\\ 0& 0& 0.01T\end{array}\right],$ 其中， $Q_0=0.001\left[\begin{array}{cc}\frac{1}{3}& \frac{1}{2}\\ \frac{1}{2}& 1\end{array}\right];$

(2)按发明内容部分步骤2所述的方法进行粒子集初始化，得到初始粒子集；

(3)按发明内容部分步骤3所述的方法获取当前时刻的量测；

(4)按发明内容部分步骤4所述的方法对粒子集进行一步预测，得到预测的粒子集；

(5)按发明内容部分步骤5所述的方法对预测的粒子集的权重进行更新；

(6)按发明内容部分步骤6所述的方法对权重更新后的粒子集进行重采样，得到当前时刻的粒子集；

(7)按发明内容部分步骤7所述的方法得到目标的状态和状态协方差估计；

(8)按发明内容部分步骤8所述的方法进行目标“状态估计-航迹”数据关联；

(9)循环执行发明内容部分步骤3～步骤8，直至雷达关机。

由附图2可以看出，在k＝13时刻，目标存在概率低于0.3，而本实施例中目标存在判决门限为γ＝0.6，基于多模粒子滤波的机动微弱目标检测前跟踪方法将出现目标漏检问题；由附图3可以看出，虽然基于多模粒子滤波的检测前跟踪方法给出了各时刻的目标状态估计，但不能给出目标航迹信息；由附图4可以看出，通过多模粒子滤波和数据关联，本发明的方法可以估计k＝13时刻漏检的目标状态，因此可以有效消除漏检问题，同时本发明方法给出了目标航迹信息。

Claims (1)

1.基于多模粒子滤波和数据关联的机动微弱目标检测前跟踪方法，其特征包括以下步骤：步骤1：变量初始化

K是总仿真时间；

T是雷达扫描周期；

N是滤波器采用的粒子数；

SNR是目标信噪比；

I_mim和I_max是目标可能的最小和最大强度；

v_min和v_max是目标可能的最小和最大速度；

γ是目标存在判决门限；

是目标初始运动模式

η₀是目标初始存在概率；

L×H表示雷达监测区域的大小；

∏_e是目标存在状态概率转移矩阵；

∏_m是目标运动模式概率转移矩阵；

Q_k是过程噪声协方差；

步骤2：令k＝0，进行粒子初始化，得到初始粒子集

其中，

包含目标的位置

速度

和强度

信息，

是目标运动模式变量，

是目标存在状态变量，i＝1，2，…，N

(1)令i＝1；

(2)目标位置初始化

$x_k^l=L\×\mathrm{rand}\left(1\right)---\left(1\right)$

$y_k^l=H\×\mathrm{rand}\left(1\right)---\left(2\right)$

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

(3)目标速度初始化

${\stackrel{\·}{x}}_k^l=v_{\min }+(v_{\max }-v_{\min })\×\mathrm{rand}\left(1\right)---\left(3\right)$

${\stackrel{\·}{y}}_k^l=v_{\min }+(v_{\max }-v_{\min })\×\mathrm{rand}\left(1\right)---\left(4\right)$

(4)目标强度初始化

$I_k^l=I_{\min }+(I_{\max }-I_{\min })\×\mathrm{rand}\left(1\right)---\left(5\right)$

(5)根据目标初始运动模式概率

初始化目标运动模式变量

(6)根据目标初始存在概率η₀初始化目标存在状态变量

(7)粒子权重初始化

$w_k^l=\frac{1}{N}---\left(6\right)$

(8)令i＝i+1，若i≤N，转(2)；

步骤3：令k＝k+1，获得k时刻的雷达量测

将雷达接收到的信号进行A/D变换，得到k时刻的雷达量测 $z_k\left\{z_k^{(p,q)}\right|p=\mathrm{1,2},...,m;q=\mathrm{1,2},...,n\},$ 送雷达数据处理计算机；

步骤4：对初始粒子集进行预测，得到粒子集其中，

(1)令i＝1；

(2)根据前一时刻的目标存在状态变量

和目标存在状态概率转移矩阵∏_e预测当前时刻目标存在状态变量

(3)若

且

利用方程(1)-方程(5)生成粒子状态并根据

随机生成目标运动模式变量

(4)若

且根据前一时刻的目标运动模式变量

和目标运动模式概率转移矩阵∏_m预测当前时刻目标运动模式变量

(5)根据

选择目标状态状态转移矩阵，其中

$F_k^{\left(1\right)}=\left[\begin{array}{ccccc}1& T& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& T& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]---\left(7\right)$

对应于目标匀速直线运动(CV模型)

$F_k^{\left(j\right)}=\left[\begin{array}{ccccc}1& \frac{\sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 0& -\frac{1-\cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 0\\ 0& \cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0& -\sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0\\ 0& \frac{1-\cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 1& \frac{\sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)}{{\mathrm{\Ω}}_k^{\left(j\right)}}& 0\\ 0& \sin \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0& \cos \left({\mathrm{\Ω}}_k^{\left(j\right)}T\right)& 0\\ 0& 0& 0& 0& 1\end{array}\right],$ j＝2，3(8)对应于目标坐标转动(CT模型)，其中

${\mathrm{\Ω}}_k^{\left(2\right)}=\frac{{-\mathrm{\α}}_m}{\sqrt{{\left({\stackrel{\·}{x}}_{k|k-1}^l\right)}^2+{\left({\stackrel{\·}{y}}_{k|k-1}^l\right)}^2}},$ ${\mathrm{\Ω}}_k^{\left(3\right)}=\frac{{\mathrm{\α}}_m}{\sqrt{{\left({\stackrel{\·}{x}}_{k|k-1}^l\right)}^2+{\left({\stackrel{\·}{y}}_{k|k-1}^l\right)}^2}}---\left(9\right)$

分别对应于顺时针坐标转动和逆时针坐标转动，α_m＞0是一个典型的机动加速度；

(6)根据目标运动模式变量

选择相应的状态转移方程，对粒子i的状态进行一步预测

$x_{k|k-1}^l=F_k^{\left(r_{k|k-1}\right)}{x}_{k-1}^l+v_k---\left(10\right)$

其中，v_k为过程噪声，其噪声协方差为Q_k；

(7)令i＝i+1，若i≤N，转(2)；

步骤5：利用当前时刻的量测对粒子权重进行更新

(1)令i＝1；

(2)确定似然区域的范围

$C_p\left(x_{k|k-1}^l\right)=\{p_0-s,...,p_0-1,p_0,p_0+1,...,p_0+s\}---\left(11\right)$

$C_q\left(y_{k|k-1}^l\right)=\{q_0-s,...,q_0-1,q_0,q_0+1,...,q_0+s\}---\left(12\right)$

其中，p₀表示最接近状态矢量元素的整数，q₀表示最接近状态矢量元素

的整数，s为预设参数；

(3)粒子权重更新

$w_{k|k-1}^l=\left\{\begin{array}{cc}\underset{p\∈C_p\left(x_{k|k-1}^l\right)}{\mathrm{\Π}}\underset{q\∈C_q\left(x_{k|k-1}^l\right)}{\mathrm{\Π}}l\left(z_k^{(p,q)}\right|x_{k|k-1}^l)& E_{k|k-1}^l=1\\ 1& E_{k|k-1}^l=0\end{array}\right.---\left(13\right)$

其中

$l\left(z_k^{(p,q)}\right|x_{k|k-1}^l)=\exp \{-\frac{h_k^{(p,q)}\left(x_{k|k-1}^l\right)(h_k^{(p,q)}\left(x_{k|k-1}^l\right)-{2z}_k^{(p,q)})}{{2\mathrm{\σ}}^2}\}---\left(14\right)$

$h_k^{(p,q)}\left(x_{k|k-1}^l\right)=\frac{{\mathrm{\Δ}}_x{\mathrm{\Δ}}_y{I}_{k|k-1}^l}{2\mathrm{\π}{\mathrm{\Σ}}^2}\exp \{-\frac{{({\mathrm{p\Δ}}_x-x_{k|k-1}^l)}^2+{({\mathrm{q\Δ}}_y-y_{k|k-1}^l)}^2}{2{\mathrm{\Σ}}^2}\}---\left(15\right)$

其中，σ²为量测噪声协方差，Δ_x和Δ_y表示传感器分辨率，∑表示传感器模糊斑点数；

(4)令i＝i+1，若i≤N，转(2)；

(5)权重归一化

$w_k^l=\frac{w_{k|k-1}^l}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{k|k-1}^j},$ i＝1，2，…，N    (16)

步骤6：重采样

对权重更新后的粒子集

进行重采样

${\{y_k^l,(1/N)\}}_{l=1}^N=\mathrm{resample}{\left\{y_{k|k-1}^l{w}_k^l\right\}}_{l=1}^N---\left(17\right)$

得到k时刻的粒子集

其中

步骤7：目标状态和状态协方差估计

(1)计算目标存在概率

${\hat{P}}_k=\frac{1}{N}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{E}_k^l---\left(16\right)$

(2)若

进行目标状态估计

${\hat{x}}_k=\frac{\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^l{E}_k^l}{\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{E}_k^l}---\left(19\right)$

和状态协方差估计

${\hat{B}}_k=\frac{1}{N}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}({\hat{x}}_k-x_k^l){({\hat{x}}_k-x_k^l)}^T---\left(20\right)$

步骤8：“状态估计-航迹”数据关联

(1)如果当前时刻不存在目标状态估计

转(4)；

(2)将目标状态估计

等效为新的量测值，状态协方差估计

等效为量测误差；

(3)若当前时刻不存在目标航迹，将新量测作为航迹头，利用m/n逻辑法进行航迹起始；否则，利用新量测对目标航迹τ_k-1或τ′_k-1进行更新，得到更新的航迹τ_k，转步骤9；

(4)如果当前时刻存在目标航迹τ_k-1或τ′_k-1，若已经连续3个时刻不存在目标状态估计，输出航迹τ_k-3并删除航迹τ′_k-1；否则，对τ_k-1或τ′_k-1进行一步预测，得到预测τ′_k；

步骤9：重复步骤3～步骤8，直至雷达关机。

Images