一种适用于闪烁噪声的多目标跟踪方法及系统与流程

本发明属于多传感器信息融合技术领域，尤其涉及一种适用于闪烁噪声的多目标跟踪方法及系统。

背景技术：

边缘分布贝叶斯滤波器是适用于杂波和噪声下多目标跟踪的有效方法，但滤波器仅适用于噪声为高斯噪声的情况，不适用于噪声为闪烁噪声的情况。如何解决闪烁噪声下多目标跟踪问题是需要探索和解决的关键技术问题。

技术实现要素：

本发明所要解决的技术问题在于提供一种适用于闪烁噪声的多目标跟踪方法与跟踪系统，旨在解决闪烁噪声下非线性运动目标的跟踪问题。

本发明是这样实现的，一种适用于闪烁噪声的多目标跟踪方法，包括：

根据前一时刻各个目标的合同分布和存在概率以及当前时刻与前一时刻的时间差，采用启发式的方法产生伽马分布的形状参数和尺度参数，进而得到当前时刻各个目标的预测合同分布和预测存在概率；

根据所述当前时刻各个目标的预测合同分布和预测存在概率，利用变分贝叶斯方法序贯处理当前时刻的测量数据，得到当前时刻各个目标的更新合同分布和更新存在概率；

利用当前时刻的测量数据生成新生目标的合同分布，并为所述新生目标指定存在概率，将所述新生目标的合同分布及存在概率分别与所述当前时刻的更新合同分布及更新存在概率进行合并，得到当前时刻各个目标的合同分布和存在概率；

从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标，并将裁减后余下目标的合同分布和存在概率作为滤波器下一次递归的输入，从所述裁减后余下的目标中提取存在概率大于第二阈值的目标，所提取出的目标的合同分布作为所述当前时刻的输出，所输出的合同分布的均值作为当前时刻目标的状态估计。

本发明还提供了一种适用于闪烁噪声的多目标跟踪系统，包括：

预测模块，用于根据前一时刻各个目标的合同分布和存在概率以及当前时刻与前一时刻的时间差，采用启发式的方法产生伽马分布的形状参数和尺度参数，进而得到当前时刻各个目标的预测合同分布和预测存在概率；

更新模块，用于根据所述当前时刻各个目标的预测合同分布和预测存在概率，利用变分贝叶斯方法序贯处理当前时刻的测量数据，得到当前时刻各个目标的更新合同分布和更新存在概率；

生成模块，用于利用当前时刻的测量数据生成新生目标的合同分布，并为所述新生目标指定存在概率，将所述新生目标的合同分布及存在概率分别与所述当前时刻的更新合同分布及更新存在概率进行合并，得到当前时刻各个目标的合同分布和存在概率；

提取模块，用于从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标，并将裁减后余下目标的合同分布和存在概率作为滤波器下一次递归的输入，从所述裁减后余下的目标中提取存在概率大于第二阈值的目标，所提取出的目标的合同分布作为所述当前时刻的输出，所输出的合同分布的均值作为当前时刻目标的状态估计。

本发明与现有技术相比，有益效果在于：本发明实施例通过用t分布建模闪烁噪声，利用变分贝叶斯方法获得闭合表达式，用各个分变量边缘分布的乘积来逼近多变量的合同分布，从而将多变量的合同估计化为对各个分变量边缘分布的迭代估计，有效解决了闪烁噪声下非线性运动目标的跟踪问题，提高多目标的跟踪精度。

附图说明

图1是本发明实施例提供的一种适用于闪烁噪声的多目标跟踪方法的流程图；

图2是本发明实施例提供的传感器70个扫描周期的测量数据；

图3是在闪烁噪声下按照本发明实施例提供的多目标跟踪方法处理得到的结果；

图4是在闪烁噪声下根据UK-PHD滤波器的多目标跟踪方法处理得到的结果；

图5是按照本发明实施例提供多目标跟踪方法与UK-PHD滤波方法在经过100次实验得到的平均OSPA距离示意图；

图6是本发明实施例提供一种适用于闪烁噪声的多目标跟踪系统的结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明实施例根据前一时刻接收的测量数据的各个目标的合同分布及存在概率预测当前时刻接收的新的测量数据的各个目标的合同分布及存在概率；根据预测的合同分布及存在概率，利用变分贝叶斯方法序贯处理当前时刻的各测量得到各个目标的更新合同分布及存在概率；分别将更新的合同分布及存在概率与新目标的合同分布及存在概率进行合并，生成当前时刻各个目标的合同分布及存在概率，这样使得本发明实施例能够很好的解决闪烁噪声下非线性运动目标的跟踪问题。

基于上述原理，本发明实施例提供了如图1所示的一种适用于闪烁噪声的多目标跟踪方法，包括：

S101，根据前一时刻各个目标的合同分布和存在概率以及当前时刻与前一时刻的时间差，采用启发式的方法产生伽马分布的形状参数和尺度参数，进而得到当前时刻各个目标的预测合同分布和预测存在概率；

S102，根据所述当前时刻各个目标的预测合同分布和预测存在概率，利用变分贝叶斯方法序贯处理当前时刻的测量数据，得到当前时刻各个目标的更新合同分布和更新存在概率；

S103，利用当前时刻的测量数据生成新生目标的合同分布，并为所述新生目标指定存在概率，将所述新生目标的合同分布及存在概率分别与所述当前时刻的更新合同分布及更新存在概率进行合并，得到当前时刻各个目标的合同分布和存在概率；

S104，从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标，并将裁减后余下目标的合同分布和存在概率作为滤波器下一次递归的输入，从所述裁减后余下的目标中提取存在概率大于第二阈值的目标，所提取出的目标的合同分布作为所述当前时刻的输出，所输出的合同分布的均值作为当前时刻目标的状态估计。

具体地，在步骤S101中，以k-1表示前一时刻，k表示当前时刻，t_k-1表示前一时刻的时间，t_k表示当前时刻的时间，所述当前时刻的测量噪声服从ξ维的t分布，以S(z_k；H_kx_k,R_k,r_k)表示所述当前时刻测量的概率密度函数，其中H_kx_k表示t分布的均值，R_k表示精度矩阵，r_k表示t分布的自由度，且前一时刻目标i的多变量合同分布为目标i的存在概率为ρ_i,k-1，其中，N表示高斯分布，g表示伽玛分布，x_i,k-1表示前一时刻第i个合同分布的状态向量，m_i,k-1表示前一时刻第i个合同分布中高斯分布的均值，P_i,k-1表示前一时刻第i个合同分布中高斯分布的方差，表示R_k的对角线元素，和γ_i,k-1表示前一时刻第i个合同分布中伽玛分布的形状参数，和η_i,k-1表示前一时刻第i个合同分布中伽玛分布的尺度参数,ξ为状态向量的维数，i＝1,2,...,N_k-1，N_k-1为前一时刻目标的总数；

根据前一时刻各个目标的合同分布和存在概率、当前时刻与前一时刻的时间差，得到当前时刻各个目标的预测合同分布当前时刻各个目标的预测存在概率为ρ_i,k|k-1＝P_s,k(t_k-t_k-1)ρ_i,k-1；其中，i＝1,2,...,N_k-1，为当前时刻第i个合同分布中高斯分布的均值，为当前时刻第i个合同分布中高斯分布的方差，Sigma点x_i,0＝m_i,k-1，Sigma点的权重l＝1,...,ξ，为目标的幸存概率，γ_i,k|k-1＝ρ_γγ_i,k-1为当前时刻第i个合同分布中伽玛分布的形状参数，η_i,k|k-1＝ρ_ηη_i,k-1为当前时刻第i个合同分布中伽玛分布的尺度参数，f为非线性函数，Q_k-1为所述接收时刻的过程噪声方差矩阵，上标T表示矩阵或向量的转置，T为采样周期，δ为已知的常数，ρ_α，ρ_β，ρ_γ，ρ_η为传播因子，取值范围为(0,1]，r_k表示自由度，为已知常数，ξ为状态向量的维数，k为一尺度参数。

在步骤S102中，设当前时刻接收到的观测集为y_k＝(y_1,k,…,y_M,k)，其中，M为当前时刻接收到测量数据总数，则所述根据所述当前时刻各个目标的预测合同分布和预测存在概率，利用变分贝叶斯方法序贯处理当前时刻的测量数据，得到当前时刻各个目标的更新合同分布和更新存在概率，包括：

S1021，以当前时刻各个目标的预测合同分布和预测存在概率作为当前时刻各个目标的初始合同分布和初始存在概率，即初始合同分布取为

初始存在概率取为其中i＝1,2,...,N_k-1，

S1022，利用变分贝叶斯方法对第1个至第M个测量数据依次进行序贯处理；

设第j个测量数据处理前各个目标的合同分布及存在概率分别为和其中，i＝1,2,...,N_k-1,1≤j≤M；由和求得用第j个测量更新时各个目标的存在概率为其中

求得用第j个测量更新时各个目标的合同分布为

其中，表示伽玛函数，tr表示矩阵的迹，表示均值向量，表示协方差矩阵，表示滤波器增益；其中Sigma点伽玛分布的形状参数为伽玛分布的尺度参数为H_k为观测矩阵，R_k为观测噪声方差矩阵，P_D,k为目标的检测概率，λ_c,k为杂波密度，I表示单位矩阵，y_j,k为当前时刻接收到的第j个测量数据，上标T表示为矩阵或向量的转置，ξ为状态向量的总维数；

若则第j个测量数据处理后目标i的合同分布为

其存在概率为其中

若则第j个测量数据处理后目标i的合同分布为

其存在概率为其中

S1023，第M个测量数据处理后的各个目标的合同分布及存在概率分别为和其中i＝1,2,...,N_k-1；将第M个测量数据处理后各个目标的合同分布及存在概率分别作为当前时刻各个目标的更新合同分布，由此得到当前时刻各个目标的更新合同分布为

及当前时刻各个目标的更新存在概率其中i＝1,…,N_k-1，

在步骤S103中，设所述当前时刻各个目标的更新合同分布为

各个目标的存在概率为ρ_i,k；其中i＝1,2,...,N_k-1，利用所述当前时刻的M个测量数据生成当前时刻新生目标的合同分布为并指定当前时刻各新生目标的存在概率为其中j＝1,2,...,M，为给定的第j个新生目标的合同分布中高斯分布的协方差，为第j个新生目标的合同分布中高斯分布的均值，由所述当前时刻第j个测量数据y_j,k＝[x_j,k y_j,k]^T产生，并且和为所述当前时刻第j个新生目标的合同分布中伽玛分布的形状参数，和为所述当前时刻第j个新生目标的合同分布中伽玛分布的尺度参数。

将所述当前时刻各个目标的更新合同分布与所述当前时刻新生目标的合同分布进行合并，得到当前时刻各个目标的合同分布为

将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并，得到所述当前时刻各个目标的存在概率为其中N_k＝N_k-1+M。

下面结合图2至图5对本实施例进行进一步地解释：

在本实施例中，考虑二维空间[-1000m,1000m]_×[-1000m,1000m]中非线性运动的目标。目标的状态由位置、速度和转弯率构成，表示为x＝[x &x y &y ω]，其中x和y分别表示位置分量，&x和&y分别表示速度分量，ω表示转弯率，上标T表示向量的转置，状态转移矩阵为过程噪声方差矩阵为△t_k＝t_k-t_k-1为当前时刻与前一时刻的时间差，σ_v和σ_w为过程噪声标准差；雷达的观测方程为：

观测噪声方差矩阵σ_r和σ_θ为观测噪声标准差，测量噪声v_k假设服从r_k＝10的t分布。

为了产生仿真数据，取幸存概率p_S,k＝1.0、检测概率p_D,k＝0.9、过程噪声的标准差σ_v＝1ms^-2，σ_w＝0.1rads^-2和观测噪声的标准差σ_r＝2m，σ_θ＝0.0003rad。一次实验的仿真观测数据如图2所示。为了处理仿真数据，将本发明实施例与闪烁噪声下的无迹卡尔曼高斯混合概率假设密度滤波器的相关参数设置为p_S,k＝1.0、p_D,k＝0.9、σ_v＝1ms^-2、σ_w＝0.1rads^-2，σ_r＝2m，σ_θ＝0.0003rad、第一阈值为10^-3、第二阈值为0.5，传播因子ρ_α＝ρ_β＝ρ_γ＝ρ_η＝0.75，伽玛分布形状参数初始值伽玛分布尺度参数初始值闪烁噪声下UK-PHD滤波器新生目标的权重w_γ＝0.1、本发明实施例的新生目标的存在概率ρ_γ＝0.1、新产生目标的协方差图3和图4为对比滤波器与本发明实施例提供的多目标跟踪方法产生的结果。将本发明实施例与现有的闪烁噪声下的UK-PHD滤波器对图2的仿真数据进行处理，100次Monte Carlo实验得到平均OSPA(Optimal Subpattern Assignment，最优亚模式分配)距离如图5所示。将现有的基于闪烁噪声的UK-PHD滤波器与本发明相比，本发明的多目标跟踪方法在闪烁噪声下对于关联不确定、检测不确定的非线性运动目标的跟踪可获得更为精确和可靠的目标状态估计、其OSPA距离比现有的这种方法得到的OSPA距离要小。

本发明还提供了如图6所示的实施例，一种多目标跟踪系统，包括：

预测模块601，用于根据前一时刻各个目标的合同分布和存在概率以及当前时刻与前一时刻的时间差，采用启发式的方法产生伽马分布的形状参数和尺度参数，进而得到当前时刻各个目标的预测合同分布和预测存在概率；

更新模块602，用于根据所述当前时刻各个目标的预测合同分布和预测存在概率，利用变分贝叶斯方法序贯处理当前时刻的测量数据，得到当前时刻各个目标的更新合同分布和更新存在概率；

生成模块603，用于利用当前时刻的测量数据生成新生目标的合同分布，并为所述新生目标指定存在概率，将所述新生目标的合同分布及存在概率分别与所述当前时刻的更新合同分布及更新存在概率进行合并，得到当前时刻各个目标的合同分布和存在概率；

提取模块604，用于从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标，并将裁减后余下目标的合同分布和存在概率作为滤波器下一次递归的输入，从所述裁减后余下的目标中提取存在概率大于第二阈值的目标，所提取出的目标的合同分布作为所述当前时刻的输出，所输出的合同分布的均值作为当前时刻目标的状态估计。

进一步地，预测模块601具体用于：

以k-1表示前一时刻，k表示当前时刻，t_k-1表示前一时刻的时间，t_k表示当前时刻的时间，所述当前时刻的测量噪声服从ξ维的t分布，以S(z_k；H_kx_k,R_k,r_k)表示所述当前时刻测量的概率密度函数，其中H_kx_k表示t分布的均值，R_k表示精度矩阵，r_k表示t分布的自由度，且前一时刻目标i的多变量合同分布为目标i的存在概率为ρ_i,k-1，其中，N表示高斯分布，g表示伽玛分布，x_i,k-1表示前一时刻第i个合同分布的状态向量，m_i,k-1表示前一时刻第i个合同分布中高斯分布的均值，P_i,k-1表示前一时刻第i个合同分布中高斯分布的方差，表示R_k的对角线元素，和γ_i,k-1表示前一时刻第i个合同分布中伽玛分布的形状参数，和η_i,k-1表示前一时刻第i个合同分布中伽玛分布的尺度参数，ξ为状态向量的维数，i＝1,2,...,N_k-1，N_k-1为前一时刻目标的总数；

根据前一时刻各个目标的合同分布和存在概率、当前时刻与前一时刻的时间差，得到当前时刻各个目标的预测合同分布当前时刻各个目标的预测存在概率为ρ_i,k|k-1＝P_s,k(t_k-t_k-1)ρ_i,k-1；其中，i＝1,2,...,N_k-1，为当前时刻第i个合同分布中高斯分布的均值，为当前时刻第i个合同分布中高斯分布的方差，Sigma点x_i,0＝m_i,k-1，Sigma点的权重l＝1,...,ξ，为目标的幸存概率，γ_i,k|k-1＝ρ_γγ_i,k-1为当前时刻第i个合同分布中伽玛分布的形状参数，η_i,k|k-1＝ρ_ηη_i,k-1为当前时刻第i个合同分布中伽玛分布的尺度参数，f为非线性函数，Q_k-1为所述接收时刻的过程噪声方差矩阵，上标T表示矩阵或向量的转置，T为采样周期，δ为已知的常数，ρ_α，ρ_β，ρ_γ，ρ_η为传播因子，取值范围为(0,1]，r_k表示自由度，为已知常数，ξ为状态向量的维数，k为一尺度参数。

进一步地，更新模块602具体用于：

以当前时刻各个目标的预测合同分布和预测存在概率作为当前时刻各个目标的初始合同分布和初始存在概率，即初始合同分布取为初始存在概率取为其中i＝1,2,...,N_k-1，

利用变分贝叶斯方法对第1个至第M个测量数据依次进行序贯处理；

设第j个测量数据处理前各个目标的合同分布及存在概率分别为和其中，i＝1,2,...,N_k-1,1≤j≤M；由和求得用第j个测量更新时各个目标的存在概率为其中

求得用第j个测量更新时各个目标的合同分布为

其中，表示伽玛函数，tr表示矩阵的迹，表示均值向量，表示协方差矩阵，表示滤波器增益；其中Sigma点伽玛分布的形状参数为伽玛分布的尺度参数为H_k为观测矩阵，R_k为观测噪声方差矩阵，P_D,k为目标的检测概率，λ_c,k为杂波密度，I表示单位矩阵，y_j,k为当前时刻接收到的第j个测量数据，上标T表示为矩阵或向量的转置，ξ为状态向量的总维数；

若则第j个测量数据处理后目标i的合同分布为

其存在概率为其中

若则第j个测量数据处理后目标i的合同分布为

其存在概率为其中

第M个测量数据处理后的各个目标的合同分布及存在概率分别为和其中i＝1,2,...,N_k-1；

将第M个测量数据处理后各个目标的合同分布及存在概率分别作为当前时刻各个目标的更新合同分布，由此得到当前时刻各个目标的更新合同分布为

及当前时刻各个目标的更新存在概率其中i＝1,…,N_k-1，

进一步地，生成模块603具体用于：

设所述当前时刻各个目标的更新合同分布为各个目标的存在概率为ρ_i,k；其中i＝1,2,...,N_k-1，利用所述当前时刻的M个测量数据生成当前时刻新生目标的合同分布为并指定当前时刻各新生目标的存在概率为其中j＝1,2,...,M，为给定的第j个新生目标的合同分布中高斯分布的协方差，为第j个新生目标的合同分布中高斯分布的均值，由所述当前时刻第j个测量数据y_j,k＝[x_j,k y_j,k]^T产生，并且和为所述当前时刻第j个新生目标的合同分布中伽玛分布的形状参数，和为所述当前时刻第j个新生目标的合同分布中伽玛分布的尺度参数。

将所述当前时刻各个目标的更新合同分布与所述当前时刻新生目标的合同分布进行合并，得到当前时刻各个目标的合同分布为

将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并，得到所述当前时刻各个目标的存在概率为其中N_k＝N_k-1+M。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。