基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法与流程
本发明属于雷达目标跟踪技术领域，特别涉及一种基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法，即基于高斯粒子势概率假设密度(cardinalized probability hypothesis density,CPHD)滤波的目标检测前跟踪方法，适用于信噪比低或数据处理量大情况下的弱小目标跟踪。
背景技术：
由于雷达反隐身和远程预警等领域的迫切需求，微弱目标探测逐渐成为目前的研究热点；当前对于微弱目标的探测使用主要使用检测前跟踪(TBD)技术，TBD技术对每次扫描不设置门限，对多帧扫描数据进行能量积累，从而对弱小目标进行检测和跟踪；TBD技术的实现方法主要有基于Hough变换方法、动态规划方法和粒子滤波方法。
随着多目标随机集研究的兴起，Mahler利用随机集PHD滤波和CPHD滤波来解决多目标跟踪问题。CPHD滤波器在传递PHD函数的同时也传递目标的势分布函数，因此比PHD滤波保留了更完整的目标数量信息，实现对目标个数的较准确估计，从而提高跟踪性能。Vo给出其序贯蒙特卡罗实现(Sequential Monte Carlo PHD，SMC-CPHD)方法。
Punithakmar K,Kirubarajan T,Sinha A.等人在其发表的A sequential Monte Carlo probability hypothesis density algorithm for multi-target track-before-detect[C].Proc.of the Signal Data Processing Small Target,2005,5913:1-8中将SMC-PHD滤波器用于红外图像多目标TBD问题中，构建了多目标状态模型和传感器量测模型，林再平，周一宇，安玮在其发表的基于势概率假设密度滤波的检测前跟踪新算法[J].红外与毫米波学报,2013,32(5):437-443中提出CPHD-TBD滤波器来增强跟踪鲁棒性。上述两种算法均采用蒙特卡洛实现，复杂环境下粒子支撑集大，导致复杂度很高。Jayesh H.K,Petar M.D.在其发表的Gaussian particle filtering[J].IEEE Trans.on Aerospace and Electronic Systems,2001:429-432中提出高斯粒子滤波器(Gaussian particle filter，GPF)，将未知状态变量的后验分布近似为高斯函数，迭代存储目标状态的均值与协方差，可以显著降低运算的复杂度，但在目标数估计上改善有限。
目前在检测和跟踪未知数目的弱小目标时面临的问题和挑战有下述问题：一方面，目标状态需要通过计算复杂的聚类方法获得，信号较弱时跟踪精度下降，另一方面，由于粒子滤波存在粒子退化、采样枯竭等问题，虽可扩大粒子支撑集以解决，但却以牺牲时间为代价，导致在实际应用中受限较大，在低信噪比下目标数的估计精度不理想。
技术实现要素：
针对上述现有技术存在的不足，本发明的目的在于提出一种基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法，该种基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法针对基于高斯粒子势概率密度假设滤波的多弱小目标检测前进行跟踪，能够大大降低运算量，实现对目标数和目标状态的准确估计，并能够有效解决解决未知数目弱小目标检测跟踪问题。
为达到上述技术目的，本发明采用如下技术方案予以实现。
一种基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法，包括以下步骤：
步骤1，初始化：令k表示k时刻，k的初始值为1，k∈{1,2,…,D}，D表示设定的最大时刻，且D为观察每个目标的运动时间；
分别将0时刻N0个目标的运动状态记为x0，将0时刻N0个目标的运动状态x0在(i,j)处的强度记为进而计算在0时刻N0个目标的运动状态x0在(i,j)处的强度条件下0时刻N0个目标的运动状态x0的后验概率密度p(x0|z0)，并分别计算0时刻N0个目标的运动状态x0的均值μ0和0时刻N0个目标的运动状态x0的协方差P0；
然后对在0时刻N0个目标的运动状态x0在(i,j)处的强度条件下0时刻N0个目标的运动状态x0的后验概率密度p(x0|z0)进行抽样，得到用于估计0时刻第p个目标的运动状态中的第ip个采样粒子表示0时刻第p个目标的运动状态中的第ip个采样粒子，Np表示0时刻第p个目标的运动状态中包含的采样粒子个数；
步骤2，确定k时刻有Nk个目标，则将k时刻第p个目标的运动状态记为其中,表示k时刻第p个目标的横坐标位置，表示k时刻第p个目标的纵坐标位置，表示k时刻第p个目标的横坐标速度，表示k时刻第p个目标的纵坐标速度，表示k时刻第p个目标的能量强度，则得到k+1时刻第p个目标的运动方程为fk(·)表示k时刻的状态转移函数，vk表示k时刻的过程处理噪声；
步骤3，确定k时刻直角坐标系下的雷达观测区域为Nx×Ny，且k时刻直角坐标系下的雷达观测区域Nx×Ny内包含Nk个目标，Nx表示k时刻直角坐标系下的雷达在x轴方向上的分辨单元数，Ny表示k时刻直角坐标系下的雷达在y轴方向上的分辨单元数，k时刻直角坐标系下的雷达观测区域中每个分辨单元(i',j')的大小为△x×△y，i'＝1,2,…,Nx,j'＝1,2,…,Ny，△x表示每个分辨单元在x轴方向上的长度，△y表示每个分辨单元在y轴方向上的长度；
根据k时刻直角坐标系下的雷达观测区域中包含Nk个目标，进而确定k时刻的红外传感器量测模型，所述k时刻的红外传感器量测模型为k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度C表示k时刻Nk个目标的扩散影响区域，所述k时刻Nk个目标的扩散影响区域为k时刻直角坐标系下的雷达观测区域Nx×Ny之外的区域；
步骤4，根据k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度计算得到在k时刻第p个目标的运动状态条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度对应的似然函数
步骤5，计算在1时刻到k-1时刻所有目标各自的量测值z1:k-1条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的后验概率密度p(xk|z1:k-1)，然后根据k+1时刻第p个目标的运动方程计算k时刻第p个目标的运动状态中第ip个采样粒子的状态预测值Np表示k时刻第p个目标的运动状态中包含的采样粒子个数；
步骤6，依次计算k-1时刻预测k时刻目标数为Nk的概率pk|k-1(Nk)和k-1时刻预测k时刻第p个目标的运动状态的强度
步骤7，计算在1时刻到k时刻所有目标各自的量测值z1:k条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的概率分布辅助函数π(xk|z1:k)，然后对所述在1时刻到k时刻所有目标各自的量测值z1:k条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的概率分布辅助函数π(xk|z1:k)进行采样并获取粒子，得到k时刻第p个目标的运动状态中第ip个采样粒子的状态估计值进而分别计算k时刻Nk个目标的状态均值估计和k时刻Nk个目标的协方差估计
步骤8，根据k时刻Nk个目标的状态均值估计和k时刻Nk个目标的协方差估计计算得到k时刻第p个的运动状态中第ip个采样粒子的权值并进行归一化操作，得到k时刻第p个的运动状态中第ip个采样粒子的归一化权值
步骤9，设置删除门限T，并将k时刻第p个的运动状态中第ip个采样粒子的归一化权值中小于设置删除门限T的高斯项删除，即删除权重小于删除门限T的高斯项，得到权重大于或等于删除门限T的k时刻第p个的运动状态中第ip个采样粒子的归一化权值同时设置合并门限U，将所述权重大于或等于删除门限T的k时刻第p个的运动状态中第ip个采样粒子的归一化权值中距离小于合并门限U的高斯项合并，得到经过删除门限T和合并门限U后k时刻第p个的运动状态中第ip个采样粒子的归一化权值
步骤10，根据k-1时刻预测k时刻目标数为Nk的概率pk|k-1(Nk)和经过删除门限T和合并门限U后k时刻第p个的运动状态中第ip个采样粒子的归一化权值计算得到k时刻目标数为Nk的概率pk(Nk)；
步骤11，根据k时刻目标数为Nk的概率pk(Nk)，计算得到k时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数Nk的估计值
步骤12，令k加1，返回步骤2，直到得到D时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数ND的估计值从而得到1时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数N1的估计值到D时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数ND的估计值
与现有技术相比，本发明具有以下优点：
第一，本发明采用高斯粒子滤波器(Gaussian particle filter，GPF)，将未知状态变量的后验分布近似为高斯函数，迭代存储目标状态的均值与协方差，可以显著降低运算的复杂度；
第二，本发明考虑检测前跟踪技术对单帧数据不进行门限判别，保留目标的原始信息，直接对原图像处理，对低信噪比下弱小目标的跟踪精度较高；
第三，本发明基于势概率密度假设(cardinalized probability hypothesis density,CPHD)滤波，降低了运算量，提高了运算效率，提高了对目标数和目标状态估计的准确性和稳定性，具有较高的时效性和工程应用性。
附图说明
下面结合附图和具体实施方式对本发明作进一步详细说明。
图1为本发明的一种基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法流程图；
图2是在信噪比为5.7dB的环境下，分别对第5、18、38、55帧记录的红外量测图；
图3(a)是在信噪比SNR＝10dB下的目标数估计示意图，水平方向为时间，单位为s，竖直方向为跟踪目标数；
图3(b)是在信噪比SNR＝7.3dB下的目标数估计示意图，水平方向为时间，单位为s，竖直方向为跟踪目标数；
图3(c)是在信噪比SNR＝5.7dB下的目标数估计示意图，水平方向为时间，单位为s，竖直方向为跟踪目标数；
图4(a)是在信噪比SNR＝10dB下的OSPA势误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA势误差距离，单位为m；
图4(b)是在信噪比SNR＝7.3dB下的OSPA势误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA势误差距离，单位为m；
图4(c)是在信噪比SNR＝5.7dB下的OSPA势误差估计，水平方向为时间，单位为s，竖直方向为OSPA势误差距离，单位为m；
图5(a)是在信噪比SNR＝10dB下的OSPA位置误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA位置误差距离，单位为m；
图5(b)是在信噪比SNR＝10dB下的OSPA位置误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA位置误差距离，单位为m；
图5(c)是在信噪比SNR＝10dB下的OSPA位置误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA位置误差距离，单位为m。
具体实施方式
参照图1，为本发明的一种基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法流程图；所述基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法，包括以下步骤：
步骤1，初始化：令k表示k时刻，k的初始值为1，k∈{1,2,…,D}，D表示设定的最大时刻，且D为观察每个目标的运动时间；本实施例中设定D＝60。
分别将0时刻N0个目标的运动状态记为x0，将0时刻N0个目标的运动状态x0在(i,j)处的强度记为并简写为z0；进而计算在0时刻N0个目标的运动状态x0在(i,j)处的强度条件下0时刻N0个目标的运动状态x0的后验概率密度p(x0|z0)，p(x0|z0)＝Ν(x0；μ0,P0)，Ν(x0；μ0,P0)表示0时刻N0个目标的运动状态x0服从均值为μ0、协方差为P0的高斯分布；N表示高斯分布函数。
其中，0时刻N0个目标的运动状态x0的均值μ0和0时刻N0个目标的运动状态x0的协方差P0，其计算表达式分别为：
其中，ip∈{1,2,...,Np}，Np表示0时刻第p个目标的运动状态包含的采样粒子个数，表示0时刻第p个目标的运动状态中第ip个采样粒子，上标T表示转置。
然后对在0时刻N0个目标的运动状态x0在(i,j)处的强度条件下0时刻N0个目标的运动状态x0的后验概率密度p(x0|z0)进行抽样，得到用于估计0时刻第p个目标的运动状态中的第ip个采样粒子表示0时刻第p个目标的运动状态中的第ip个采样粒子，Np表示0时刻第p个目标的运动状态中包含的采样粒子个数。
其中，估计0时刻第p个目标的运动状态需要Np个采样粒子，1≤Np≤+∞，Np逼近任何形式的概率密度分布，且Np的取值越趋于正无穷，估计的目标运动状态越接近真实情况；考虑运算效率和复杂度等因素，本实施例采样粒子数选取为Np＝2000。
步骤2，建立目标运动模型：确定k时刻有Nk个目标，则将k时刻第p个目标的运动状态记为其中,表示k时刻第p个目标的横坐标位置，表示k时刻第p个目标的纵坐标位置，表示k时刻第p个目标的横坐标速度，表示k时刻第p个目标的纵坐标速度，表示k时刻第p个目标的能量强度，则得到k+1时刻第p个目标的运动方程为fk(·)表示k时刻的状态转移函数，vk表示k时刻的过程处理噪声。
步骤3，建立红外传感器量测模型：确定k时刻直角坐标系下的雷达观测区域为Nx×Ny，且k时刻直角坐标系下的雷达观测区域Nx×Ny内包含Nk个目标，Nx表示k时刻直角坐标系下的雷达在x轴方向上的分辨单元数，Ny表示k时刻直角坐标系下的雷达在y轴方向上的分辨单元数，k时刻直角坐标系下的雷达观测区域中每个分辨单元(i',j')的大小为△x×△y，i'＝1,2,…,Nx,j'＝1,2,…,Ny，△x表示每个分辨单元在x轴方向上的长度，△y表示每个分辨单元在y轴方向上的长度。
根据k时刻直角坐标系下的雷达观测区域中包含Nk个目标，进而确定k时刻的红外传感器量测模型，所述k时刻的红外传感器量测模型为k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度其表达式为：
其中，nk(i,j)表示k时刻Nk个目标在分辨单元(i,j)处的量测噪声，且所述量测噪声为独立高斯白噪声；C表示k时刻Nk个目标的扩散影响区域，所述k时刻Nk个目标的扩散影响区域为k时刻直角坐标系下的雷达观测区域Nx×Ny之外的区域；表示k时刻第p个目标的运动状态在分辨单元(i,j)处的信号强度扩展函数，其表达式近似为：
其中，Σ表示点扩散方差，△x表示每个分辨单元在x轴方向上的长度，△y表示每个分辨单元在y轴方向上的长度，表示k时刻第p个目标的能量强度，表示k时刻第p个目标的横坐标位置，表示k时刻第p个目标的纵坐标位置。
步骤4，根据k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度计算得到在k时刻第p个目标的运动状态条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度对应的似然函数其表达式为：
其中，表示k时刻第p个目标的运动状态在x轴方向上的扩散区域范围，表示k时刻第p个目标的运动状态在y轴方向上的扩散区域范围，q表示k时刻第p个目标的扩散影响区域半径，(r,s)表示k时刻第p个目标的运动状态中扩散影响区域最大的分辨单元，r表示k时刻第p个目标的运动状态中扩散影响区域最大的分辨单元在x轴方向上的分辨单元数，s表示k时刻第p个目标的运动状态中扩散影响区域最大的分辨单元在y轴方向上的分辨单元数；表示在k时刻第p个目标的运动状态条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度对应的第p个目标扩散影响区域内的似然函数，表示在k时刻第p个目标的运动状态条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度对应的第p个目标扩散影响区域外的似然函数；∈表示属于，表示不属于；且将k时刻第p个目标的信噪比记为其表达式为：Σ表示点扩散方差，lg表示对数操作，σ表示k时刻Nk个目标在分辨单元(i,j)处的量测噪声nk(i,j)的方差；本实施例中将红外图像中k时刻第p个目标的信噪比低于10dB以下时视为弱目标，并将扩散影响区域小于6×6的对应目标视为小目标。
步骤5，计算在1时刻到k-1时刻所有目标各自的量测值z1:k-1条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的后验概率密度p(xk|z1:k-1)，其表达式为：p(xk|z1:k-1)＝Ν(xk-1；μk-1,Pk-1)，Ν(xk-1；μk-1,Pk-1)表示k-1时刻直角坐标系下雷达观测区域中包含的Nk-1个目标的运动状态xk-1均值为μk-1、协方差为Pk-1的高斯分布；所述在1时刻到k-1时刻所有目标各自的量测值条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的后验概率密度p(xk|z1:k-1)为根据1时刻到k-1时刻所有目标各自的量测值z1:k-1得到k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的概率；其中，当k＝1时z1:k-1表示k-1时刻到1时刻所有目标各自的量测值。
然后根据k+1时刻第p个目标的运动方程计算k时刻第p个目标的运动状态中第ip个采样粒子的状态预测值其表达式为：fk(·)表示k时刻的状态转移函数，xk-1表示k-1时刻Nk-1个目标的运动状态，μk-1表示k-1时刻Nk-1个目标的运动状态xk-1服从高斯分布时的均值，Pk-1表示k-1时刻Nk-1个目标的运动状态xk-1服从高斯分布时的协方差，表示k时刻第p个目标的运动状态中第ip个采样粒子，表示k+1时刻第p个目标的运动方程，Np表示k时刻第p个目标的运动状态中包含的采样粒子个数。
步骤6，依次计算k-1时刻预测k时刻目标数为Nk的概率pk|k-1(Nk)和k-1时刻预测k时刻第p个目标的运动状态的强度
所述k-1时刻预测k时刻目标数为Nk的概率pk|k-1(Nk)，其表达式为：
其中，Nk表示k时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数，j表示k时刻直角坐标系下的雷达观测区域Nx×Ny内存活的目标个数，m-j表示k时刻直角坐标系下的雷达观测区域Nx×Ny内的新生目标公式，l表示k-1时刻直角坐标系下的雷达观测区域Nx×Ny内的存活目标个数，l-j表示k时刻直角坐标系下的雷达观测区域Nx×Ny内消亡目标的个数；所述存活目标为直角坐标系下的雷达观测区域Nx×Ny内k-1时刻存在且k时刻也存在的目标，所述新生目标为直角坐标系下的雷达观测区域Nx×Ny内k-1时刻不存在但k时刻存在的目标，所述消亡目标为直角坐标系下的雷达观测区域Nx×Ny内k-1时刻存在但k时刻不存在的目标。
pΓ,k(·)表示k时刻新生目标的概率，下标Γ表示新生目标，下标s表示存活目标，ps,k(·)表示k时刻存活目标的概率，→表示趋近于，l表示第l个目标，ps,k-1(l)表示k-1时刻第l个目标在直角坐标系下的雷达观测区域Nx×Ny内的存活功率，表示k时刻第j个目标在直角坐标系下的雷达观测区域Nx×Ny内的存活功率，ps,k表示k时刻任意一个目标的存活概率。
步骤7，根据k-1时刻预测k时刻目标数为Nk的概率pk|k-1(Nk)，计算k-1时刻预测k时刻第p个目标的运动状态的强度其表达式为：
其中，表示k-1时刻预测k时刻仍然存活的第p个目标的强度，即不会随着时间消亡；表示k时刻第p个新生目标的运动状态的强度，ps,k表示k时刻任意一个目标的存活概率，表示k-1时刻第p个目标的运动状态中第ip'个采样粒子的权值，表示k时刻第p个目标的运动状态，表示k-1时刻预测k时刻第p个存活目标的运动状态中第ip个采样粒子的高斯项权值，表示k-1时刻预测k时刻第p个存活目标的运动状态中第ip个采样粒子的高斯项协方差，表示k时刻第p个新生目标的运动状态中第ip”个采样粒子的高斯项权值，表示k时刻第p个新生目标的运动状态中第ip”个采样粒子的高斯项均值，表示k时刻第p个新生目标的运动状态中第ip”个采样粒子的高斯项协方差，JΓ,k表示k时刻每一个新生目标的运动状态中包含的采样粒子个数，JΓ,k-1表示k-1时刻每一个新生目标的运动状态中包含的采样粒子个数。
步骤8，计算在1时刻到k时刻所有目标各自的量测值z1:k条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的概率分布辅助函数π(xk|z1:k)，π(xk|z1:k)＝p(xk|z1:k-1)，p(xk|z1:k-1)表示在1时刻到k-1时刻所有目标各自的量测值z1:k-1条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的后验概率密度，z1:k-1表示1时刻到k-1时刻所有目标各自的量测值，xk表示k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态。
然后对所述在1时刻到k时刻所有目标各自的量测值z1:k条件下k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk的概率分布辅助函数π(xk|z1:k)进行采样并获取粒子，得到k时刻第p个目标的运动状态中第ip个采样粒子的状态估计值进而分别计算k时刻Nk个目标的状态均值估计和k时刻Nk个目标的协方差估计其表达式分别为：
其中，表示k时刻第p个目标的运动状态中第ip'个采样粒子的权值，上标T表示转置。
步骤9，根据k时刻Nk个目标的状态均值估计和k时刻Nk个目标的协方差估计计算得到k时刻第p个的运动状态中第ip个采样粒子的权值并进行归一化操作，得到k时刻第p个的运动状态中第ip个采样粒子的归一化权值其表达式分别为：
其中，Ν(xk；μk,Pk)表示k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk均值为μk、协方差为Pk的高斯分布，表示k时刻第p个目标的运动状态中第ip个采样粒子条件下zk的后验概率密度，表示在1时刻到k时刻所有目标各自的量测值z1:k条件下k时刻第p个目标的运动状态中第ip个采样粒子的概率分布辅助函数，表示k时刻第p个目标的运动状态中第ip个采样粒子，zk表示k时刻直角坐标系下雷达观测区域中包含的Nk个目标的运动状态xk在(i,j)处的强度
步骤10，高斯项随着时间的递增而增加，为了避免其数目无限增长，设置删除门限T，并将k时刻第p个的运动状态中第ip个采样粒子的归一化权值中小于设置删除门限T的高斯项删除，即删除权重小于删除门限T的高斯项，得到权重大于或等于删除门限T的k时刻第p个的运动状态中第ip个采样粒子的归一化权值同时设置合并门限U，将所述权重大于或等于删除门限T的k时刻第p个的运动状态中第ip个采样粒子的归一化权值中距离小于合并门限U的高斯项合并，得到经过删除门限T和合并门限U后k时刻第p个的运动状态中第ip个采样粒子的归一化权值进而分别计算k时刻Nk个目标的状态均值优化估计和k时刻Nk个目标的协方差优化估计并将所述k时刻Nk个目标的状态均值优化估计和k时刻Nk个目标的协方差优化估计记为k时刻直角坐标系下雷达观测区域中包含的Nk个目标各自的目标状态，所述目标状态包含目标的位置和速度；本实施例中T＝10-5，U＝4。
步骤11，根据k-1时刻预测k时刻目标数为Nk的概率pk|k-1(Nk)和经过删除门限T和合并门限U后k时刻第p个的运动状态中第ip个采样粒子的归一化权值计算得到k时刻目标数为Nk的概率pk(Nk)，其表达式为：
其中，g(·)表示分布概率转移函数，表示经过删除门限T和合并门限U后k时刻第p个的运动状态中第ip个采样粒子的归一化权值，z1:k表示在1时刻到k时刻所有目标各自的量测值。
步骤12，根据k时刻目标数为Nk的概率pk(Nk)，计算得到k时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数Nk的估计值其表达式为：
其中，argmax表示反向取最大值运算。
步骤13，令k加1，返回步骤2，直到得到D时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数ND的估计值从而得到1时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数N1的估计值到D时刻直角坐标系下的雷达观测区域Nx×Ny内包含的目标个数ND的估计值
通过以下仿真实验对本发明效果作进一步验证说明。
(一)仿真实验条件：
仿真在MATLAB7.0软件下进行。
仿真场景中设置目标采用匀速运动模型，且k+1时刻Nk个目标的运动方程为Xk+1＝FXk+Rk，其中Rk为零均值高斯白噪声，F表示状态转移函数；目标检测概率为1，生存概率为0.99，目标出生概率为0.1，目标消失概率为0.01，目标衍生概率为0。整个检测与跟踪过程持续时间为60帧，帧间间隔为1秒，像素分辨单元△x＝△y＝1，图像大小为80×80的序列图像，传感器模糊系数∑＝0.7，目标强度服从Ι∈U(18,22)的均匀分布。修剪门限T为10-5，合并门限U为4。
设置三个目标的初始状态分别为：
X1＝[5,0.15,6.4,0.1,20]T,X2＝[7.8,0.12,1,0.1,20]T,X3＝[2,0.1,8.2,0.12,20]T。目标1存在于整个60帧监测时间区间，目标2在t＝10s进入量测区域，t＝40s消失，目标3在t＝20s进入量测区域，t＝60s消失。
(二)仿真实验内容：
目标信噪比和目标检测与跟踪性能息息相关，故仿真在不同信噪比环境下进行实验，信噪比环境分别为10dB，7.3dB，5.7dB。每帧新生粒子数2000，并进行100次蒙特卡洛仿真实验，然后分别使用基于粒子势概率假设密度滤波的目标检测前跟踪方法(SMC-CPHD-TBD)、基于高斯粒子概率假设密度滤波的目标检测前跟踪方法(GPF-PHD-TBD)和本发明的基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法(GPF-CPHD-TBD)三种算法进行比较。
仿真1：仿真场景中设置在信噪比为5.7dB的环境下，分别对第5、18、38、55帧记录的红外量测图像，像素分辨单元△x＝△y＝1，图像大小为80×80的序列图像，传感器模糊系数∑＝0.7；参照图2，为在信噪比为5.7dB的环境下，分别对第5、18、38、55帧记录的红外量测图。
仿真2:在未知数目的弱小目标跟踪过程的，对目标数目的估计是体现跟踪算法精度的关键。仿真2设置在信噪比SNR分别在10dB、7.3dB和5.7dB的情况下，对三种跟踪方法(SMC-CPHD-TBD，GPF-CPHD-TBD和GPF-PHD-TBD)进行目标数目跟踪比较，其中，GPF-CPHD-TBD方法是本发明所提方法。如图3(a)和图3(b)所示，图3(a)是在信噪比SNR＝10dB下的目标数估计示意图，水平方向为时间，单位为s，竖直方向为跟踪目标数；图3(b)是在信噪比SNR＝7.3dB下的目标数估计示意图，水平方向为时间，单位为s，竖直方向为跟踪目标数。比较可知，在SNR＝10dB和7.3dB时，使用三种方法均可较为准确地对弱小目标进行目标数量估计，本发明方法比其余两种方法更为精准，如图3(c)所示，为信噪比SNR＝5.7dB下的目标数估计示意图，水平方向为时间，单位为s，竖直方向为跟踪目标数。在较差环境(SNR＝5.7dB)中，SMC-CPHD-TBD方法和GPF-PHD-TBD方法跟踪性能明显下降，而本发明所提方法依旧能准确估计目标数。
仿真3：在考虑估计多个弱小目标状态的同时，也需要估计目标数目，选用最优子模式分配(OSPA)距离作为评价准则指标，OSPA包括位置和集合势两部分距离，其中距离敏感性参数p表征距离误差，水平调节数c表征集合势误差，考虑距离误差比势误差更为重要。设定OSPA参数c＝5，p＝2，并假定新生目标位置未知，在每帧的观测图像中要寻找新生目标。仿真3针对OSPA势误差进行试验，OSPA势误差是描绘方法跟踪目标数目和原本真实目标数目之间的误差值，OSPA值越小，表明跟踪方法的性能越好。图4(a)是在信噪比SNR＝10dB下的OSPA势误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA势误差距离，单位为m；图4(b)是在信噪比SNR＝7.3dB下的OSPA势误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA势误差距离，单位为m；图4(c)是在信噪比SNR＝5.7dB下的OSPA势误差估计，水平方向为时间，单位为s，竖直方向为OSPA势误差距离，单位为m。通过比较三种算法的OSPA势误差，体现了在目标数估计方面的总体性能。从结果图分析可知，随着信噪比的不断降低，三种方法的OSPA值均有所升高，尤其在SNR＝5.7dB时增幅明显，意味着三种方法估计性能随信噪比降低都有不同程度的下降。但相较于另两种方法，本发明所提方法GPF-CPHD-TBD的OSPA势误差值始终较低，即使在信噪比为5.7dB的较差环境中跟踪性能仍保持良好，表明所提方法可以稳定捕捉目标，准确估计目标数目。
仿真4：仿真4针对OSPA距离误差进行试验，OSPA距离误差是描绘算法跟踪目标的坐标位置信息和原本真实目标位置坐标信息之间的误差值，OSPA值越小，表明跟踪方法的性能越好；参照图5(a)至图5(c)，图5(a)是在信噪比SNR＝10dB下的OSPA位置误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA位置误差距离，单位为m；图5(b)是在信噪比SNR＝10dB下的OSPA位置误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA位置误差距离，单位为m；图5(c)是在信噪比SNR＝10dB下的OSPA位置误差估计示意图，水平方向为时间，单位为s，竖直方向为OSPA位置误差距离，单位为m；通过比较三种方法OSPA位置误差，体现目标跟踪精度的总体性能。SMC-CPHD-TBD方法随着目标数目增加、信噪比下降等相关情况，跟踪性能下降较为明显，这是因为传统的粒子滤波方法本身存在粒子退化等现象，而重采样由可能带来采样枯竭等问题。分析比较GPF-PHD-TBD方法和本发明所提GPF-CPHD-TBD方法，前者OSPA位置误差值随信噪比的下降升高明显，表明对目标位置的估计精准性下降幅度较大。而所提方法的OSPA位置误差值在三种信噪比环境中均保持最低，特别是低信噪比下仍性能良好，说明比其余两种算法更精准估计目标状态，保持更稳定的跟踪性能，鲁棒性更佳。
仿真5：进行三种方法的时间复杂度实验，表1给出了三种方法的时间复杂度的比较分析，即SMC-CPHD-TBD方法、GPF-CPHD-TBD方法和GPF-PHD-TBD方法。
表1
本发明所提GPF-CPHD-TBD方法约为SMC-CPHD-TBD方法所耗时的38％，因为GPF采用拟蒙特卡洛高斯粒子采样，从而可以得到更为均匀的粒子数，每个高斯项只许分配少量粒子就可实现目标较为准确的跟踪。虽然GPF-CPHD-TBD方法比GPF-PHD-TBD方法耗时稍长，但前者比后者更能稳定地估计目标数目，从而实现对目标发现和状态的准确估计，这对低信噪比条件下的弱小目标跟踪研究意义重大。
综合上述处理结果，可以发现本发明方法用高斯粒子滤波的思想，迭代更新目标状态的均值与协方差，无需重采样，并结合CPHD的TBD算法实现对多弱小目标进行检测和跟踪；仿真实验表明，使用本发明方法有效地降低了运算量，提高了运算效率，提高了对目标数和目标状态估计的准确性和稳定性，解决了在低信噪比下未知数目的弱小目标检测跟踪问题，具有较高的时效性和工程应用性。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
- 还没有人留言评论。精彩留言会获得点赞！