CN104766320A - 阈值化量测下的多伯努利滤波弱目标检测与跟踪方法 - Google Patents

Abstract

本发明公开了一种阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，主要解决现有技术在低信噪比情况下弱小目标跟踪效率低、跟踪性能不理想、计算复杂度高的问题。本发明实现的方法是：先设定一个量测阈值，提取灰度图像数据中超过量测阈值的像素点的索引，将索引集作为量测集进行存储。再利用设定的量测阈值，计算灰度图像中各个像素点的检测概率。根据存储的量测集和像素点的检测概率，计算量测似然。采用多伯努利滤波对目标进行检测与跟踪。本发明在低信噪比情况下，能有效地提高弱小目标的跟踪效率和跟踪精度，降低计算复杂度和存储空间。

Description

阈值化量测下的多伯努利滤波弱目标检测与跟踪方法

技术领域

本发明属于制导技术领域，更进一步涉及目标跟踪技术领域中的一种阈值化量测下的多伯努利滤波弱目标检测与跟踪方法。本发明通过计算量测阈值化后的似然函数，采用高斯粒子多伯努利检测前跟踪方法，实现对未知数目的多个弱小目标的精确检测与跟踪。

背景技术

弱小目标的检测与跟踪是红外预警系统，精确制导系统，卫星遥感系统中的一项关键技术，在现代化高科技发展中具有重要的应用价值。传统的弱小目标检测算法采用的是先检测后跟踪(DBT)方法，它需要对每帧数据设置固定阈值或自适应阈值，此类算法结构简单，易于实现，在信噪比较高时有较好的效果。然而在低信噪比下，很难找到一个合适的阈值来分割目标。近年来出现的一种检测前跟踪(TBD)方法，很好地解决了低信噪比下的弱目标检测跟踪问题。它不对数据做单帧阈值检测，而是将未经阈值处理的数据全部输入跟踪方法进行跟踪，待跟踪若干帧之后，目标能量积累到一定程度再做目标有无的判断。

哈尔滨工程大学在其拥有的专利“基于粒子滤波的红外弱小目标检测与跟踪方法及装置”(专利申请号：2012101631401，授权公布号CN 102722706B)中公开了一种基于粒子滤波的检测前跟踪方法来检测和跟踪弱小目标。首先，对原始图像进行背景抑制，得到含有目标和噪声的去背景图像；其次，根据设定的阈值和初始目标出现概率,在全视场范围内采样粒子进行目标状态的预测和更新，在状态更新过程中计算粒子权重并进行归一化；再次，根据归一化后的粒子权重对粒子进行重采样；最后，进行目标检测，得到目标状态的估计值。该方法存在的不足之处在于，传统的粒子滤波采样粒子容量大以及重采样带来的粒子枯竭，导致在低信噪比情况下目标跟踪性能不理想的问题。

Cuiyun Li,Hongbing Ji,Qibing Zou,Sujun Wang在其发表的论文“A NOVELMULTI-BERNOULLI FILTER FOR JOINT TARGET DETECTION AND TRACKING”(Information and Communication and Communications Technologies(IETICT2013)，IET International Conference on 2013,Pages：176-180)中公开了一种高斯粒子多伯努利检测前跟踪方法。该方法将基于随机有限集的高斯粒子MeMBer滤波用于多个弱小目标检测与跟踪中，它是粒子滤波多伯努利检测前跟踪方法的一种推广。首先，将目标随机集用一个多伯努利随机集来近似；其次，进行多伯努利状态预测和状态更新，此过程不需要进行粒子重采样，只需少量的采样粒子计算目标的状态均值和方差，降低了存储空间，提高了跟踪精度；再次，进行伯努利项的修剪与合并；最后，进行目标状态估计值的提取。该方法的不足之处在于，传统的弱目标检测方法需要存储和传输每帧图像的全部数据，存在计算量大，计算复杂度高以及高存储空间占用问题，并且雷达，声呐，被动红外传感器大多是阈值传感器，使得该方法在实际应用中还存在缺陷。

发明内容

本发明的目的在于针对上述现有技术的不足，提出了一种阈值化量测下的多伯努利滤波弱目标检测与跟踪方法。

实现本发明的基本思路是：首先，利用虚警概率，对灰度图像数据设定量测阈值；其次，利用设置的量测阈值，提取每帧灰度图像数据中超过阈值的像素点的索引，并将提取的索引集作为量测集进行存储；再次，利用设置的量测阈值，计算每帧灰度图像中各个像素点的检测概率；最后，根据存储的量测集和像素点的检测概率，计算量测似然，使用多伯努利滤波对目标进行检测与跟踪。

实现本发明的具体步骤如下：

(1)初始化目标状态：

(1a)根据目标运动场景，设定初始时刻存在目标的状态参数，包括目标个数、位置、速度、强度、存在概率以及目标的运动状态方差；

(1b)用上述设定的参数作为目标的初始分布，并用高斯粒子多伯努利随机集的参数集形式表示；

(2)预测目标状态：

(2a)利用前一时刻目标多伯努利随机集的后验分布和状态方程，预测当前时刻存活目标的多伯努利随机集；

(2b)根据目标运动场景，设定新生目标的个数、位置、速度、强度、存在概率以及目标的运动状态方差；

(2c)用上述设定的参数作为当前时刻新生目标的预测分布，并用高斯粒子多伯努利随机集的参数集形式表示；

(3)更新目标状态：

(3a)根据像素点的虚警概率设定量测阈值，利用当前时刻的灰度图像数据，提取所有强度超过量测阈值的像素点的索引，并将索引集作为当前时刻的量测集；

(3b)根据量测方程，计算当前时刻灰度图像中每个像素点强度的概率分布，利用设定的量测阈值，计算每个像素点强度超过量测阈值的概率，作为当前时刻灰度图像中像素点的检测概率；

(3c)利用当前时刻的量测集和当前时刻灰度图像中每个像素点的检测概率，计算当前时刻目标随机集的量测似然值；

(3d)利用当前时刻目标随机集的量测似然值，更新当前时刻预测完成的目标多伯努利随机集，得到当前时刻目标多伯努利随机集的后验分布；

(4)修剪与合并：

(4a)剔除当前时刻更新后的目标多伯努利随机集中存在概率低于0.001的目标；

(4b)合并当前时刻更新后的目标多伯努利随机集中，欧氏距离小于2倍的目标影响邻域参数的目标，合并之后在(0，1)之间重新设置目标的存在概率；

(5)状态提取：

(5a)利用当前时刻修剪与合并之后的目标多伯努利随机集，统计所有存在概率大于0.5的目标数目，将统计结果作为当前时刻的真实目标数；

(5b)利用当前时刻修剪与合并之后的目标多伯努利随机集，提取所有存在概率大于0.5的目标状态均值，将提取结果作为当前时刻真实存在的目标状态；

(6)判断所有灰度图像是否处理完毕，若是，执行步骤(7)，否则，执行步骤(2)，处理下一帧图像；

(7)结束。

本发明与现有技术相比，具有以下优点：

第一，由于本发明采用了对灰度图像数据设定量测阈值，提取所有强度超过阈值的像素点的索引，并将索引集作为当前时刻的量测集，克服了现有技术中计算量大，计算复杂度高以及高存储空间占用问题，使得本发明具有低信噪比情况下弱小目标跟踪效率高的优点。

第二，由于本发明采用了高斯粒子多伯努利滤波方法，用高斯分布来表示目标的概率密度，用量测阈值化后的似然函数更新目标多伯努利随机集，克服了现有技术中粒子滤波带来的低信噪比情况下弱目标跟踪性能不理想的问题，使得本发明检测和跟踪多个弱小目标的精度得到提高。

第三，由于本发明采用了对量测阈值化后的数据进行检测前跟踪，满足了雷达，声呐，被动红外传感器大多是阈值传感器的要求，克服了现有技术不易被实际应用的问题，使得本发明可被应用于现实的目标跟踪技术中。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图对本发明作进一步的详细说明。

参照附图1，本发明的具体步骤如下。

步骤1，初始化目标状态。

根据目标运动场景，设定初始时刻存在目标的状态参数，包括目标个数、位置、速度、强度、存在概率以及目标的运动状态方差。用上述设定的参数作为目标的初始分布，并用高斯粒子多伯努利随机集的参数集形式表示。

令初始时刻k＝0，目标初始分布用高斯粒子多伯努利随机集的参数集表示形式如下：

${\left\{(r_0^{\left(i\right)},m_0^{\left(i\right)},P_0^{\left(i\right)})\right\}}_{i=1}^{M_0}$

其中，表示初始时刻第i个目标的存在概率，表示一个用5维向量表示的目标运动状态均值，其中x₀表示初始时刻目标的横坐标，表示初始时刻目标的水平速度，y₀表示初始时刻目标的纵坐标，表示初始时刻目标的垂直速度，In₀表示初始时刻目标强度，T表示转置操作，表示初始时刻目标状态方差，M₀表初始时刻的目标数目。

步骤2，预测目标状态。

目标在k时刻的状态可以用一个5维向量表示其中(x_k,y_k)、和In_k分别表示k时刻目标的位置、速度和强度。当目标距离观测站非常远时，可近似认为其做匀速运动，运动方程为：

x_k＝Fx_k-1+w_k

其中，F表示状态转移矩阵，w_k表示状态噪声。

假定k-1(k＞＝1)时刻的目标多伯努利随机集后验分布为：

${\left\{(r_{k-1}^{\left(i\right)},m_{k-1}^{\left(i\right)},P_{k-1}^{\left(i\right)})\right\}}_{i=1}^{M_{k-1}}$

新生目标的多伯努利随机集为：

${\left\{(r_{B,k}^{\left(i\right)},m_{B,k}^{\left(i\right)},P_{B,k}^{\left(i\right)})\right\}}_{i=1}^{M_{B,k}}$

则k时刻预测的目标多伯努利随机集为：

${\left\{(r_{A,k|k-1}^{\left(i\right)},m_{A,k|k-1}^{\left(i\right)},P_{A,k|k-1}^{\left(i\right)})\right\}}_{i=1}^{M_{k-1}}\∪{\left\{(r_{B,k}^{\left(i\right)},m_{B,k}^{\left(i\right)},P_{B,k}^{\left(i\right)})\right\}}_{i=1}^{M_{B,k}}$

其中，表示k-1时刻第i个目标的存在概率，表示k-1时刻第i个目标的状态均值，表示k-1时刻第i个目标的状态方差，M_k-1表示k-1时刻的目标数目；Β表示新生目标，表示k时刻第i个新生目标的预测存在概率，表示k时刻第i个新生目标的预测状态均值，表示k时刻第i个新生目标的预测状态方差，M_Β,k表示k时刻新生目标的目标数目；Α表示存活目标，表示从k-1时刻到k时刻第i个存活目标的预测存在概率，表示从k-1时刻到k时刻第i个存活目标的预测状态均值，表示从k-1时刻到k时刻第i个存活目标的预测状态方差。

利用k-1时刻目标多伯努利随机集的后验分布和状态方程，预测k时刻存活目标的多伯努利随机集。

具体的预测方法可由下列步骤来完成。

第1步，对k-1时刻目标多伯努利随机集后验分布，按照下式，进行高斯采样抽取高斯采样粒子：

$x_{k-1}^{\left(i\right)\left(j\right)}~N(m_{k-1}^{\left(i\right)},P_{k-1}^{\left(i\right)})$

其中，表示k-1时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子的状态，～表示高斯采样操作，表示k-1时刻第i个目标的后验概率密度，表示k-1时刻第i个目标的状态均值，表示k-1时刻第i个目标的状态方差。

第2步，将高斯采样粒子代入下述的状态方程中，预测每个高斯采样粒子的状态值：

$x_k^{\left(i\right)\left(j\right)}=F{x}_{k-1}^{\left(i\right)\left(j\right)}+w_k$

其中，表示k时刻预测后的第i个目标的第j个高斯采样粒子的状态，F表示状态转移矩阵，表示k-1时刻第i个目标的第j个高斯采样粒子的状态，w_k表示k时刻的状态噪声。

第3步，按照下式，预测存活目标的状态均值：

$m_{A,k|k-1}^{\left(i\right)}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^{\left(i\right)\left(j\right)}$

其中，表示从k-1时刻到k时刻第i个存活目标的预测状态均值，Α表示存活目标，N表示高斯采样粒子数目，表示k时刻预测后的第i个目标的第j个高斯采样粒子的状态。

第4步，按照下式，预测存活目标的状态方差：

$P_{A,k|k-1}^{\left(i\right)}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}[x_k^{\left(i\right)\left(j\right)}-m_{A,k|k-1}^{\left(i\right)}]{[x_k^{\left(i\right)\left(j\right)}-m_{A,k|k-1}^{\left(i\right)}]}^T$

其中，表示从k-1时刻到k时刻第i个存活目标的预测状态方差，Α表示存活目标，N表示高斯采样粒子数目，表示k时刻预测后的第i个目标的第j个高斯采样粒子的状态，表示从k-1时刻到k时刻第i个存活目标的预测状态均值，T表示转置操作。

第5步，按照下式，预测存活目标的存在概率：

$r_{P,k|k-1}^{\left(i\right)}=p_s{r}_{k-1}^{\left(i\right)}$

其中，表示从k-1时刻到k时刻第i个存活目标的预测存在概率，Α表示存活目标，p_S表示目标存活的概率，表示k-1时刻第i个目标的存在概率。

根据目标运动场景，设定新生目标的个数、位置、速度、强度、存在概率以及目标的运动状态方差；用上述设定的参数作为k时刻新生目标的预测分布，并用高斯粒子多伯努利随机集的参数集形式表示。

k时刻新生目标预测分布用高斯粒子多伯努利随机集的参数集形式表示如下：

${\left\{(r_{B,k}^{\left(i\right)},m_{B,k}^{\left(i\right)},P_{B,k}^{\left(i\right)})\right\}}_{i=1}^{M_{B,k}}$

其中，表示k时刻第i个新生目标的预测存在概率，Β表示新生目标，表示一个用5维向量表示的k时刻第i个新生目标的预测状态均值，其中x_k表示k时刻预测新生目标的横坐标，表示k时刻预测新生目标的水平速度，y_k表示k时刻预测新生目标的纵坐标，表示k时刻预测新生目标的垂直速度，In_k表示k时刻预测新生目标强度，T表示转置操作，表示k时刻第i个新生目标预测状态方差，M_Β,k表示k时刻预测新生目标数目。

步骤3，更新目标状态。

考虑一个凝视型红外成像系统，观测模型就是k时刻给定目标状态x_k下在像素点l处的强度y_l：

$y_l=\left\{\begin{array}{cc}{h}_l\left(x_k\right)+v_l& l\∈T\left(x_k\right)\\ v_l& l\∉T\left(x_k\right)\end{array}\right.$

其中，y_l表示k时刻灰度图像中像素点l的强度，x_k表示k时刻的目标状态，h_l(x_k)表示k时刻目标状态x_k对像素点l处量测的贡献，v_l表示量测噪声，T(x_k)表示受目标x_k影响到的邻域像素的坐标集合。

根据像素点的虚警概率设定量测阈值，利用k时刻的灰度图像数据，提取所有强度超过量测阈值的像素点的索引，并将索引集作为k时刻的量测集。

按照下式，根据虚警概率，设定量测阈值：

$\mathrm{\β}={\∫}_{\mathrm{\γ}}^{\∞}N(0,{\mathrm{\σ}}^2)\mathrm{dy}$

其中，β表示像素点的虚警概率，γ表示设定的量测阈值，y表示像素点的强度，σ²表示量测噪声方差，N(0,σ²)表示均值为0，方差为σ²的高斯分布。

按照下式，利用设定的量测阈值，计算k时刻的量测集：

Z_k＝{l|y_l＞γ,l＝1,...,imgN}

其中，Z_k表示k时刻的量测集，y_l表示像素点l的强度，γ表示设定的量测阈值，imgN表示k时刻灰度图像中像素点的个数。

根据量测方程，计算k时刻灰度图像中每个像素点强度的概率分布，利用设定的量测阈值，计算每个像素点强度超过量测阈值的概率，作为k时刻灰度图像中像素点的检测概率。

根据量测方程，计算k时刻给定目标状态x下像素点l处的概率分布：

$p\left(y_l\right|x_k)=\left\{\begin{array}{cc}N(h_l\left(x_k\right),{\mathrm{\σ}}^2)& l\∈T\left(x_k\right)\\ N(0,{\mathrm{\σ}}^2)& l\∉T\left(x_k\right)\end{array}\right.$

其中，p(y_l|x_k)表示k时刻灰度图像中像素点l的强度的概率分布，y_l表示k时刻灰度图像中像素点l的强度，x_k表示k时刻的目标状态，h_l(x_k)表示k时刻目标状态x_k对像素点l处量测的贡献，T(x_k)表示受目标x_k影响到的邻域像素坐标集合，σ²表示高斯分布的方差，N(h_l(x_k),σ²)表示均值为h_l(x_k)，方差为σ²的高斯分布，N(0,σ²)表示均值为0，方差为σ²的高斯分布。

按照下式，计算k时刻灰度图像中像素点的检测概率：

$P_l={\∫}_{\mathrm{\γ}}^{\∞}p\left(y_l\right|x_k){\mathrm{dy}}_l$

其中，P_l表示k时刻灰度图像中像素点l的检测概率，y_l表示当前时刻灰度图像中像素点l的强度，p(y_l|x_k)表示像素点l的强度的概率分布，γ表示设定的量测阈值。

对于受目标影响的像素点l∈T(x_k)，按照下式，求得像素点l的检测概率α_l：

${\mathrm{\α}}_l={\∫}_{\mathrm{\γ}}^{\∞}p\left(y_l\right|x_k){\mathrm{dy}}_l={\∫}_{\mathrm{\γ}}^{\∞}N(h_l\left(x_k\right),{\mathrm{\σ}}^2){\mathrm{dy}}_l$

对于不受目标影响的像素点按照下式，求得像素点l的虚警概率β_l：

${\mathrm{\β}}_l={\∫}_{\mathrm{\γ}}^{\∞}p\left(y_l\right|x_k){\mathrm{dy}}_l={\∫}_{\mathrm{\γ}}^{\∞}N(0,{\mathrm{\σ}}^2){\mathrm{dy}}_l$

利用k时刻的量测集和k时刻灰度图像中每个像素点的检测概率，计算k时刻目标随机集的量测似然值。

根据上述模型假设，k时刻目标随机集的量测似然函数为：

$f\left(Z_k\right|X)=f\left(Z_k\right)\underset{x_k\∈X}{\mathrm{\Π}}{g}_{Z_k}\left(x_k\right)$

其中， $f\left(Z_k\right)=\underset{l}{\overset{\mathrm{imgN}}{\mathrm{\Π}}}(1-{\mathrm{\β}}_l)\underset{l\∈Z_k}{\mathrm{\Π}}\frac{{\mathrm{\β}}_l}{1-{\mathrm{\β}}_l},g_{Z_k}\left(x_k\right)=\underset{l\∈\mathrm{\Ω}}{\mathrm{\Π}}\frac{1-{\mathrm{\α}}_l}{1-{\mathrm{\β}}_l}\underset{l\∈\mathrm{\Ψ}}{\mathrm{\Π}}\frac{{\mathrm{\α}}_l}{{\mathrm{\β}}_l},$ X表示k时刻的目标随机集，Z_k表示k时刻的量测集，f(Z_k|X)表示k时刻目标随机集X的量测似然，x_k表示k时刻目标随机集X中的目标状态，表示目标状态x_k的似然比函数，Ω＝T(x_k)\Z_k，Ω表示一差集，\表示求差集操作，T(x_k)表示受目标x_k影响到的邻域像素坐标集合，Ψ＝T(x_k)∩Z_k，Ψ表示一并集，∩表示求并集操作，α_l表示像素点l受目标x_k影响时的检测概率，β_l表示像素点l不受目标x_k影响时的虚警概率，imgN表示像素点的个数。

假设k时刻预测的目标多伯努利随机集为：

${\left\{(r_{k|k-1}^{\left(i\right)},m_{k|k-1}^{\left(i\right)},P_{k|k-1}^{\left(i\right)})\right\}}_{i=1}^{M_{k|k-1}}$

则更新的目标多伯努利随机集后验分布为：

${\left\{(r_k^{\left(i\right)},m_k^{\left(i\right)},P_k^{\left(i\right)})\right\}}_{i=1}^{M_k}$

其中，表示从k-1时刻到k时刻预测完成后的第i个目标的存在概率，表示从k-1时刻到k时刻预测完成后的第i个目标的状态均值，表示从k-1时刻到k时刻预测完成后的第i个目标的状态方差，M_k|k-1表示k-1时刻到k时刻预测完成的目标数目；表示k时刻更新完成后的第i个目标的存在概率，表示k时刻更新完成后的第i个目标的状态均值，表示k时刻更新完成后的第i个目标的状态方差，M_k表示k时刻更新完成后的目标数目。

利用k时刻目标随机集的量测似然值，更新k时刻预测完成的目标多伯努利随机集，得到k时刻目标多伯努利随机集的后验分布。

具体的更新方法可由下列步骤来完成。

第1步，按照下式，对k时刻预测的目标状态，进行高斯采样抽取高斯采样粒子：

$x_k^{\left(i\right)\left(j\right)}~N(m_{k|k-1}^{\left(i\right)},P_{k|k-1}^{\left(i\right)})$

其中，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，～表示高斯采样操作，表示从k-1时刻到k时刻预测完成后的第i个目标的状态均值，表示从k-1时刻到k时刻预测完成后的第i个目标的状态方差，表示从k-1时刻到k时刻预测完成后的第i个目标的概率密度。

第2步，按照下式，计算每个高斯采样粒子的权值：

$w_k^{\left(i\right)\left(j\right)}=\frac{1}{N}{g}_{Z_k}\left(x_k^{\left(i\right)\left(j\right)}\right)$

其中，表示k时刻第i个目标的第j个高斯采样粒子的权值，N表示高斯采样粒子数目，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，Z_k表示k时刻的量测集，表示采样粒子状态的似然比函数。

第3步，按照下式，更新k时刻目标状态均值：

$m_k^{\left(i\right)}=\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}{x}_k^{\left(i\right)\left(j\right)}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}$

其中，表示k时刻第i个目标的更新状态均值，表示k时刻第i个目标的第j个高斯采样粒子的权值，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，N表示高斯采样粒子数目。

第4步，按照下式，更新k时刻目标状态方差：

$P_k^{\left(i\right)}=\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}[x_k^{\left(i\right)\left(j\right)}-m_k^{\left(i\right)}]{[x_k^{\left(i\right)\left(j\right)}-m_k^{\left(i\right)}]}^T}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}$

其中，表示k时刻第i个目标的更新状态方差，表示k时刻第i个目标的第j个高斯采样粒子的权值，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，表示k时刻第i个目标的更新状态均值，N表示高斯采样粒子数目，T表示转置操作。

第5步，按照下式，更新k时刻目标存在概率：

$r_k^{\left(i\right)}=\frac{r_{k|k-1}^{\left(i\right)}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}{1-r_{k|k-1}^{\left(i\right)}+r_{k|k-1}^{\left(i\right)}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}$

其中，表示k时刻第i个目标的更新存在概率，表示从k-1时刻到k时刻预测完成后的第i个目标的存在概率，N表示高斯采样粒子数目，表示k时刻第i个目标的第j个高斯采样粒子的权值。

步骤4，修剪与合并。

剔除k时刻更新后的目标多伯努利随机集中存在概率低于0.001的目标；合并目标多伯努利随机集中欧氏距离小于2倍的目标影响邻域参数的目标，合并之后在(0，1)之间重新设置目标的存在概率。

步骤5，状态提取。

利用k时刻修剪与合并之后的目标多伯努利随机集，统计所有存在概率大于0.5的目标数目，将统计结果作为k时刻的真实目标数；提取所有存在概率大于0.5的目标状态均值，将提取结果作为k时刻真实存在的目标状态。

步骤6，判断所有灰度图像是否处理完毕，若是，执行步骤7，否则，令k＝k+1，返回执行步骤2，处理下一帧图像。

步骤7，结束。

下面结合附图2的仿真，对本发明的效果做进一步说明。

1，仿真条件：

本发明在Intel(R)Core(TM)i5-4430CPU3.00GHz 3.00GHz处理器的电脑上，采用MATLAB R2014a软件完成仿真。

仿真场景设置：为了验证本发明提出的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法能精确的检测和跟踪弱小目标，本发明的仿真实验场景为在的二维空间[0,80]×[0,80]内有4个做匀速直线运动的目标，相继出现和消失在不同的时间帧，整个时间持续40帧，目标1存在时间是t＝1s到t＝34s，目标2存在时间是从t＝10s到t＝38s，目标3存在时间是从t＝5s到t＝40s，目标4存在时间是从t＝15s到t＝37s。附图2(a)为目标真实轨迹，其中横轴表示目标的横坐标，纵轴表示目标的纵坐标，轨迹中的圆圈o表示目标的起始点。

目标状态方程为：

x_k＝Fx_k-1+w_k

其中， $F=\left[\begin{array}{ccc}{F}_x& & \\ & F_x& \\ & & 1\end{array}\right],F_x=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right],$ w_k表示协方差矩阵为Q的零均值高斯白噪声 $Q=\left[\begin{array}{ccc}{Q}_s& & \\ & Q_s& \\ & & q_{2}T\end{array}\right],Q_s=\left[\begin{array}{cc}\frac{q_1}{3}{T}^3& \frac{q_1}{2}{T}^2\\ \frac{q_1}{2}{T}^2& q_{2}T\end{array}\right].$ T表示采样周期，q₁，q₂分别表示目标运动噪声和目标强度噪声的功率谱密度。

量测方程为：

$y_l=\left\{\begin{array}{cc}{h}_l\left(x_k\right)+v_l& l\∈T\left(x_k\right)\\ v_l& l\∉T\left(x_k\right)\end{array}\right.$

其中， $h_l\left(x_k\right)=\frac{{\mathrm{\Δ}}_x{\mathrm{\Δ}}_y{\mathrm{In}}_k}{2\mathrm{\π}{\mathrm{\σ}}_h^2}\exp (-\frac{{({\mathrm{\Δ}}_{x}i-p_{x,k})}^2+{({\mathrm{\Δ}}_{y}j-p_{y,k})}^2}{2{\mathrm{\σ}}_h^2}),$ v_l表示量测噪声，服从均值为0，方差为σ²的高斯分布，σ_h表示传感器的模糊系数，(p_x,k,p_y,k)表示目标x_k的位置，i，j分别表示像素点l的横纵坐标，In_k表示k时刻目标的强度。

多伯努利滤波器初始化时给定第一个目标的准确位置m₀＝[5,1,70,-1,In]^T，新生目标的位置给定准确值，即r_Γ＝0.02， $m_{\mathrm{\Γ}}^{\left(1\right)}={[\mathrm{25,1,37,1},\mathrm{In}]}^T,m_{\mathrm{\Γ}}^{\left(2\right)}={[\mathrm{40,1,45},-1,\mathrm{In}]}^T,m_{\mathrm{\Γ}}^{\left(3\right)}={[\mathrm{40,1,28,1},\mathrm{In}]}^T,$ P_Γ＝diag([0.5,0.5,0.5,0.5,0.2]²)。

仿真参数设置：Δ_x＝Δ_y＝1，采样周期T＝1，目标运动噪声功率谱密度q₁＝0.001，目标强度噪声功率谱密度q₂＝0.01，传感器的模糊系数σ_h＝1，目标状态强度In＝30，每个高斯项采样的粒子数目L＝100，存在概率低于0.001的伯努利项被去除，最大伯努利项数目T_max＝20，目标的存活概率p_s＝0.99，量测阈值取为经验值γ＝2，目标影响的邻域参数为2，则伯努利项的合并门限为4。

根据信噪比计算公式：

$\mathrm{SNR}=10\log {\left[\frac{\mathrm{In}{\mathrm{\Δ}}_x{\mathrm{\Δ}}_y}{2\mathrm{\π}{\mathrm{\σ}}_h^2\mathrm{\σ}}\right]}^2$

当信噪比SNR＝9时，由上式求得量测噪声标准差σ＝1.7，再根据像素点的虚警概率公式：

$\mathrm{\β}={\∫}_{\mathrm{\γ}}^{\∞}N(0,{\mathrm{\σ}}^2)\mathrm{dy}$

当设置像素点的虚警概率β＝0.1时，由上式求得量测阈值γ≈2.2。

为了证明仿真效果，在信噪比SNR＝9时进行50次的蒙特卡罗仿真，计算估计的目标数目均值和脱靶距离OSPA，OSPA的参数为c＝10，p＝2。附图2(b)是目标跟踪轨迹，图2(b)中的横轴表示应用本发明方法跟踪结果的横坐标，纵轴表示应用本发明方法跟踪结果的纵坐标，轨迹中的黑点表示应用本发明方法的跟踪结果。附图2(c)是平均估计目标数目，图2(c)中的横轴表示时间，纵轴表示目标数目，图2(c)中实线表示目标真实数目，黑点表示应用本发明方法平均估计目标数目，虚线表示目标数目的估计方差。附图2(d)是目标位置跟踪精度，图2(d)中的横轴表示时间，纵轴表示OSPA位置误差。附图2(e)是目标数目跟踪精度，图2(e)中的横轴表示时间，纵轴表示OSPA势误差。附图2(f)是目标整体跟踪精度，图2(f)中的横轴表示时间，纵轴表示整体OSPA误差。

2，仿真内容与结果分析：

附图2(b)给出了应用本发明实现的目标跟踪轨迹图，由图2(b)可知，本发明的方法在信噪比较低时完成了对未知数目的多个弱小目标的精确检测和跟踪。附图2(c)给出了目标数的估计均值，由图2(c)可知，本发明的方法准确的估计了场景中每一时刻存在的目标数，并未发生多估或者少估的状况。附图2(d)给出了目标跟踪结果的位置精度，并用脱靶距离OSPA位置误差进行表示，由图2(d)可知，OSPA位置误差比较小并且整体比较稳定。图2(e)给出了目标跟踪结果的目标数目精度，并用脱靶距离OSPA势误差进行表示。由图2(e)可知，当有目标新生或者目标消失时，对应时间点为t＝5s，10s，15s，35s，40s，由于不能立刻检测到目标的状态变化，导致OSPA势误差出现较大起伏，但是当目标全部出现并稳定存在时，OSPA势误差趋于稳定。附图2(f)给出了目标跟踪结果的整体OSPA误差，由图2(f)可知，整体OSPA误差比较小，其波动起伏高的时刻为t＝5s，10s，15s，35s，40s，主要原因是这些时刻存在目标新生或者目标的消失，其变化趋势与OSPA势误差变化趋势相一致。

综上所述，从仿真效果图的分析可知，本发明提出的一种阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，实现了对低信噪比下未知数目的多个弱小目标的精确检测和跟踪。目标数估计准确，目标跟踪精度高，跟踪性能良好，同时由于采用了量测阈值化，使得量测集存储的是像素点的索引，大大降低了计算的复杂度和储存空间，并且针对雷达，声呐，被动红外等传感器大多是阈值传感器的应用局限，该发明在现实工程应用中更有优势。

Claims (8)

1.一种阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，包括如下步骤：

(1)初始化目标状态：

(1a)根据目标运动场景，设定初始时刻存在目标的状态参数，包括目标个数、位置、速度、强度、存在概率以及目标的运动状态方差；

(1b)用上述设定的参数作为目标的初始分布，并用高斯粒子多伯努利随机集的参数集形式表示；

(2)预测目标状态：

(2a)利用前一时刻目标多伯努利随机集的后验分布和状态方程，预测当前时刻存活目标的多伯努利随机集；

(2b)根据目标运动场景，设定新生目标的个数、位置、速度、强度、存在概率以及目标的运动状态方差；

(2c)用上述设定的参数作为当前时刻新生目标的预测分布，并用高斯粒子多伯努利随机集的参数集形式表示；

(3)更新目标状态：

(3a)根据像素点的虚警概率设定量测阈值，利用当前时刻的灰度图像数据，提取所有强度超过量测阈值的像素点的索引，并将索引集作为当前时刻的量测集；

(3b)根据量测方程，计算当前时刻灰度图像中每个像素点强度的概率分布，利用设定的量测阈值，计算每个像素点强度超过量测阈值的概率，作为当前时刻灰度图像中像素点的检测概率；

(3c)利用当前时刻的量测集和当前时刻灰度图像中每个像素点的检测概率，计算当前时刻目标随机集的量测似然值；

(3d)利用当前时刻目标随机集的量测似然值，更新当前时刻预测完成的目标多伯努利随机集，得到当前时刻目标多伯努利随机集的后验分布；

(4)修剪与合并：

(4a)剔除当前时刻更新后的目标多伯努利随机集中存在概率低于0.001的目标；

(4b)合并当前时刻更新后的目标多伯努利随机集中，欧氏距离小于2倍的目标影响邻域参数的目标，合并之后在(0，1)之间重新设置目标的存在概率；

(5)状态提取：

(5a)利用当前时刻修剪与合并之后的目标多伯努利随机集，统计所有存在概率大于0.5的目标数目，将统计结果作为当前时刻的真实目标数；

(5b)利用当前时刻修剪与合并之后的目标多伯努利随机集，提取所有存在概率大于0.5的目标状态均值，将提取结果作为当前时刻真实存在的目标状态；

(6)判断所有灰度图像是否处理完毕，若是，执行步骤(7)，否则，执行步骤(2)，处理下一帧图像；

(7)结束。

2.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(1b)中所述的高斯粒子多伯努利随机集的参数集表示形式如下：

${\left\{(r_0^{\left(i\right)},m_0^{\left(i\right)},P_0^{\left(i\right)})\right\}}_{i=1}^{M_0}$

其中，表示初始时刻第i个目标的存在概率，表示一个用5维向量表示的目标运动状态均值，其中x₀表示初始时刻目标的横坐标，表示初始时刻目标的水平速度，y₀表示初始时刻目标的纵坐标，表示初始时刻目标的垂直速度，In₀表示初始时刻目标强度，T表示转置操作，表示初始时刻目标状态方差，M₀表初始时刻的目标数目。

3.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(2a)中所述的当前时刻存活目标的多伯努利随机集的预测方法如下：

第1步，按照下式，对k-1时刻目标多伯努利随机集后验分布，进行高斯采样抽取高斯采样粒子：

$X_{k-1}^{\left(i\right)\left(j\right)}~N(m_{k-1}^{\left(i\right)},P_{k-1}^{\left(i\right)})$

其中，表示k-1时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子的状态，～表示高斯采样操作，表示k-1时刻第i个目标的后验概率密度，表示k-1时刻第i个目标的状态均值，表示k-1时刻第i个目标的状态方差；

第2步，将高斯采样粒子代入下述的状态方程中，预测每个高斯采样粒子的状态值：

$X_k^{\left(i\right)\left(j\right)}={\mathrm{Fx}}_{k-1}^{\left(i\right)\left(j\right)}+w_k$

其中，表示k时刻预测后的第i个目标的第j个高斯采样粒子的状态，F表示状态转移矩阵，表示k-1时刻第i个目标的第j个高斯采样粒子的状态，w_k表示k时刻的状态噪声；

第3步，按照下式，预测存活目标的状态均值：

$m_{A,k|k-1}^{\left(i\right)}=\frac{1}{N}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{x}_k^{\left(i\right)\left(j\right)}$

其中，表示从k-1时刻到k时刻第i个存活目标的预测状态均值，Α表示存活目标，N表示高斯采样粒子数目，表示k时刻预测后的第i个目标的第j个高斯采样粒子的状态；

第4步，按照下式，预测存活目标的状态方差：

$P_{A,k|k-1}^{\left(i\right)}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}[x_k^{\left(i\right)\left(j\right)}-m_{A,k|k-1}^{\left(i\right)}]{[x_k^{\left(i\right)\left(j\right)}-m_{A,k|k-1}^{\left(i\right)}]}^T$

其中，表示从k-1时刻到k时刻第i个存活目标的预测状态方差，Α表示存活目标，N表示高斯采样粒子数目，表示k时刻预测后的第i个目标的第j个高斯采样粒子的状态，表示从k-1时刻到k时刻第i个存活目标的预测状态均值，T表示转置操作；

第5步，按照下式，预测存活目标的存在概率：

$r_{A,k|k-1}^{\left(i\right)}=p_s{r}_{k-1}^{\left(i\right)}$

其中，表示从k-1时刻到k时刻第i个存活目标的预测存在概率，Α表示存活目标，p_S表示目标存活的概率，表示k-1时刻第i个目标的存在概率。

4.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(2c)中所述的用高斯粒子多伯努利随机集的参数集表示形式如下：

${\left\{(r_{B,k}^{\left(i\right)},m_{B,k}^{\left(i\right)},P_{B,k}^{\left(i\right)})\right\}}_{i=1}^{M_{B,k}}$

其中，表示k时刻第i个新生目标的预测存在概率，Β表示新生目标，表示一个用5维向量表示的k时刻第i个新生目标的预测状态均值，其中x_k表示k时刻预测新生目标的横坐标，表示k时刻预测新生目标的水平速度，y_k表示k时刻预测新生目标的纵坐标，表示k时刻预测新生目标的垂直速度，In_k表示k时刻预测新生目标强度，T表示转置操作，表示k时刻第i个新生目标预测状态方差，M_Β,k表示k时刻预测新生目标数目。

5.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(3a)中所述计算当前时刻的量测集方法如下：

第1步，按照下式，计算设定的量测阈值：

$\mathrm{\β}={\∫}_{\mathrm{\γ}}^{\∞}N(0,{\mathrm{\σ}}^2)\mathrm{dy}$

其中，β表示像素点的虚警概率，γ表示设定的量测阈值，y表示像素点的强度，σ²表示量测噪声方差，N(0,σ²)表示均值为0，方差为σ²的高斯分布；

第2步，按照下式，计算k时刻的量测集：

Z_k＝{l|y_l＞γ,l＝1,...,imgN}

其中，Z_k表示k时刻的量测集，y_l表示像素点l的强度，γ表示设定的量测阈值，imgN表示k时刻灰度图像中像素点的个数。

6.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(3b)中所述计算当前时刻灰度图像中像素点的检测概率方法如下：

第1步，按照下式，建立量测方程：

$y_l=\left\{\begin{array}{cc}{h}_l\left(x_k\right)+v_l& l\∈T\left(x_k\right)\\ v_l& l\∉T\left(x_k\right)\end{array}\right.$

其中，y_l表示k时刻灰度图像中像素点l的强度，x_k表示k时刻的目标状态，h_l(x_k)表示k时刻目标状态x_k对像素点l处量测的贡献，v_l表示量测噪声，T(x_k)表示受目标x_k影响到的邻域像素的坐标集合；

第2步，按照下式，计算当前时刻灰度图像中像素点强度的概率分布：

$p\left(y_l\right|x_k)=\left\{\begin{array}{cc}N(h_l\left(x_k\right),{\mathrm{\σ}}^2)& l\∈T\left(x_k\right)\\ N(0,{\mathrm{\σ}}^2)& l\∉T\left(x_k\right)\end{array}\right.$

其中，p(y_l|x_k)表示k时刻灰度图像中像素点l的强度的概率分布，y_l表示k时刻灰度图像中像素点l的强度，x_k表示k时刻的目标状态，h_l(x_k)表示k时刻目标状态x_k对像素点l处量测的贡献，T(x_k)表示受目标x_k影响到的邻域像素坐标集合，σ²表示高斯分布的方差，N(h_l(x_k),σ²)表示均值为h_l(x_k)，方差为σ²的高斯分布，N(0,σ²)表示均值为0，方差为σ²的高斯分布；

第3步，按照下式，计算当前时刻灰度图像中像素点的检测概率：

$P_l={\∫}_{\mathrm{\γ}}^{\∞}p\left(y_l\right|x_k){\mathrm{dy}}_l$

其中，P_l表示k时刻灰度图像中像素点l的检测概率，y_l表示当前时刻灰度图像中像素点l的强度，p(y_l|x_k)表示像素点l的强度的概率分布，γ表示设定的量测阈值。

7.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(3c)中所述计算当前时刻目标随机集的量测似然公式如下：

$f\left(Z_k\right|X)=f\left(Z_k\right)\underset{x_k\∈X}{\mathrm{\Π}}{g}_{Z_k}\left(x_k\right)$

其中， $f\left(Z_k\right)=\underset{l}{\overset{\mathrm{imgN}}{\mathrm{\Π}}}(1-{\mathrm{\β}}_l)\underset{l\∈Z_k}{\mathrm{\Π}}\frac{{\mathrm{\β}}_l}{1-{\mathrm{\β}}_l},g_{Z_k}\left(x_k\right)=\underset{l\∈\mathrm{\Ω}}{\mathrm{\Π}}\frac{1-{\mathrm{\α}}_l}{1-{\mathrm{\β}}_l}\underset{l\∈\mathrm{\Ψ}}{\mathrm{\Π}}\frac{{\mathrm{\α}}_l}{{\mathrm{\β}}_l},$ X表示k时刻的目标随机集，Z_k表示k时刻的量测集，f(Z_k|X)表示k时刻目标随机集X的量测似然，x_k表示k时刻目标随机集X中的目标状态，表示目标状态x_k的似然比函数，Ω＝T(x_k)\Z_k，Ω表示一差集，\表示求差集操作，T(x_k)表示受目标x_k影响到的邻域像素坐标集合，Ψ＝T(x_k)∩Z_k，Ψ表示一并集，∩表示求并集操作，α_l表示像素点l受目标x_k影响时的检测概率，β_l表示像素点l不受目标x_k影响时的虚警概率，imgN表示像素点的个数。

8.根据权利要求1所述的阈值化量测下的多伯努利滤波弱目标检测与跟踪方法，其特征在于，步骤(3d)中所述计算当前时刻目标多伯努利随机集的后验分布方法如下：

第1步，按照下式，对k时刻预测的目标状态，进行高斯采样抽取高斯采样粒子：

$x_k^{\left(i\right)\left(j\right)}~N(m_{k|k-1}^{\left(i\right)},P_{k|k-1}^{\left(i\right)})$

其中，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，～表示高斯采样操作，表示从k-1时刻到k时刻预测完成后的第i个目标的状态均值，表示从k-1时刻到k时刻预测完成后的第i个目标的状态方差，表示从k-1时刻到k时刻预测完成后的第i个目标的概率密度；

第2步，按照下式，计算每个高斯采样粒子的权值：

$w_k^{\left(i\right)\left(j\right)}=\frac{1}{N}{g}_{Z_k}\left(x_k^{\left(i\right)\left(j\right)}\right)$

其中，表示k时刻第i个目标的第j个高斯采样粒子的权值，N表示高斯采样粒子数目，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，Z_k表示k时刻的量测集，表示采样粒子状态的似然比函数；

第3步，按照下式，更新k时刻目标状态均值：

$m_k^{\left(i\right)}=\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}{x}_k^{\left(i\right)\left(j\right)}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}$

其中，表示k时刻第i个目标的更新状态均值，表示k时刻第i个目标的第j个高斯采样粒子的权值，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，N表示高斯采样粒子数目；

第4步，按照下式，更新k时刻目标状态方差：

$P_k^{\left(i\right)}=\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}[x_k^{\left(i\right)\left(j\right)}-m_k^{\left(i\right)}]{[x_k^{\left(i\right)\left(j\right)}-m_k^{\left(i\right)}]}^T}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}$

其中，表示k时刻第i个目标的更新状态方差，表示k时刻第i个目标的第j个高斯采样粒子的权值，表示k时刻对第i个目标进行高斯采样抽取的第j个高斯采样粒子状态，表示k时刻第i个目标的更新状态均值，N表示高斯采样粒子数目，T表示转置操作；

第5步，按照下式，更新k时刻目标存在概率：

$r_k^{\left(i\right)}=\frac{r_{k|k-1}^{\left(i\right)}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}{1-r_{k|k-1}^{\left(i\right)}+r_{k|k-1}^{\left(i\right)}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)\left(j\right)}}$

其中，表示k时刻第i个目标的更新存在概率，表示从k-1时刻到k时刻预测完成后的第i个目标的存在概率，N表示高斯采样粒子数目，表示k时刻第i个目标的第j个高斯采样粒子的权值。