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

摘要

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

说明书

技术领域

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

背景技术

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

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

Cuiyun Li,Hongbing Ji,Qibing Zou,Sujun Wang在其发表的论文“A NOVEL  MULTI-BERNOULLI FILTER FOR JOINT TARGET DETECTION AND TRACKING” (Information and Communication and Communications Technologies(IETICT 2013)，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}}\cup {\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}{\Sigma }}{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}{\Sigma }}[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\in T\left(x_k\right)\\ v_l& l\notin 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时刻的量测集。

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

$\beta ={\int }_{\gamma }^{\infty }N(0,{\sigma }^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),{\sigma }^2)& l\in T\left(x_k\right)\\ N(0,{\sigma }^2)& l\notin 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={\int }_{\gamma }^{\infty }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：

${\alpha }_l={\int }_{\gamma }^{\infty }p\left(y_l\right|x_k){\mathrm{dy}}_l={\int }_{\gamma }^{\infty }N(h_l\left(x_k\right),{\sigma }^2){\mathrm{dy}}_l$

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

${\beta }_l={\int }_{\gamma }^{\infty }p\left(y_l\right|x_k){\mathrm{dy}}_l={\int }_{\gamma }^{\infty }N(0,{\sigma }^2){\mathrm{dy}}_l$

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

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

$f\left(Z_k\right|X)=f\left(Z_k\right)\underset{x_k\in X}{\Pi }{g}_{Z_k}\left(x_k\right)$

其中，$f\left(Z_k\right)=\underset{l}{\overset{\mathrm{imgN}}{\Pi }}(1-{\beta }_l)\underset{l\in Z_k}{\Pi }\frac{{\beta }_l}{1-{\beta }_l},g_{Z_k}\left(x_k\right)=\underset{l\in \Omega }{\Pi }\frac{1-{\alpha }_l}{1-{\beta }_l}\underset{l\in \Psi }{\Pi }\frac{{\alpha }_l}{{\beta }_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}{\Sigma }}{w}_k^{\left(i\right)\left(j\right)}{x}_k^{\left(i\right)\left(j\right)}}{\underset{j=1}{\overset{N}{\Sigma }}{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}{\Sigma }}{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}{\Sigma }}{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}{\Sigma }}{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}{\Sigma }}{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-4430CPU@3.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\in T\left(x_k\right)\\ v_l& l\notin T\left(x_k\right)\end{array}\right)$

其中，$h_l\left(x_k\right)=\frac{{\Delta }_x{\Delta }_y{\mathrm{In}}_k}{2\pi {\sigma }_h^2}\exp (-\frac{{({\Delta }_{x}i-p_{x,k})}^2+{({\Delta }_{y}j-p_{y,k})}^2}{2{\sigma }_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_{\Gamma }^{\left(1\right)}={[\mathrm{25,1,37,1},\mathrm{In}]}^T,m_{\Gamma }^{\left(2\right)}={[\mathrm{40,1,45},-1,\mathrm{In}]}^T,m_{\Gamma }^{\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}{\Delta }_x{\Delta }_y}{2\pi {\sigma }_h^2\sigma }\right]}^2$

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

$\beta ={\int }_{\gamma }^{\infty }N(0,{\sigma }^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势误差变化趋势相一致。

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