CN106022340A - 一种改进的高斯混合势概率假设密度滤波方法 - Google Patents
一种改进的高斯混合势概率假设密度滤波方法 Download PDFInfo
- Publication number
- CN106022340A CN106022340A CN201610325272.8A CN201610325272A CN106022340A CN 106022340 A CN106022340 A CN 106022340A CN 201610325272 A CN201610325272 A CN 201610325272A CN 106022340 A CN106022340 A CN 106022340A
- Authority
- CN
- China
- Prior art keywords
- target
- hypothesis density
- probability hypothesis
- probability
- potential distribution
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 26
- 238000001914 filtration Methods 0.000 title claims abstract description 19
- 238000009826 distribution Methods 0.000 claims abstract description 88
- 238000005259 measurement Methods 0.000 claims abstract description 51
- 239000011159 matrix material Substances 0.000 claims abstract description 33
- 238000006243 chemical reaction Methods 0.000 claims abstract description 9
- 238000011156 evaluation Methods 0.000 claims abstract description 8
- 239000000203 mixture Substances 0.000 claims description 17
- 230000004083 survival effect Effects 0.000 claims description 8
- 239000013598 vector Substances 0.000 claims description 6
- 238000001514 detection method Methods 0.000 claims description 5
- 238000012546 transfer Methods 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 abstract description 6
- 238000009966 trimming Methods 0.000 abstract 1
- 238000004422 calculation algorithm Methods 0.000 description 15
- 101001086191 Borrelia burgdorferi Outer surface protein A Proteins 0.000 description 6
- 238000004088 simulation Methods 0.000 description 5
- 230000009466 transformation Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 206010027476 Metastases Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 230000008034 disappearance Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000009401 metastasis Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000013138 pruning Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/30—Noise filtering
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Operations Research (AREA)
- Multimedia (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种改进的高斯混合势概率假设密度滤波方法，步骤如下：1，组成目标的状态集和目标强度函数；2，初始化初始目标的概率假设密度及势分布；3，对目标状态集在k+1时刻的概率假设密度及势分布进行预测，得到k+1时刻的概率假设密度及势分布；4，对目标状态集在k+1时刻的概率假设密度及势分布进行更新，得到k+1时刻的概率假设密度及势分布，对真实协方差矩阵和真实偏差进行无偏转换，设置椭球门限对量测集合进行简化，减小当前观测集中的观测数目；5，对目标的强度函数的高斯项进行修剪合并，提取目标状态估计进行性能评估；6，重复步骤3～5，跟踪目标直到目标消失。本发明有利于雷达数据信息的直接运用，且减小了滤波器的计算量。
Description
一、技术领域
本发明涉及目标跟踪技术领域，特别是一种改进的高斯混合势概率假设密度滤波方法。
二、背景技术
随着非传统防务及安全挑战的不断涌现，多目标跟踪算法的研究成为一个热点。而在多目标跟踪算法中，有两类主要方法，数据关联(如PDA、JPDA)和绕过关联直接处理(如PHD、CPHD)。在目标个数较多且含杂波时，数据关联(PDA等)的计算量会非常大，不利于工程应用。而近年来，多目标跟踪研究专家Ronald P.S.Mahler教授提出了基于有限集统计学(FISST)的随机有限集(random finite set，RFS)理论，在此基础上推出概率假设密度(Probability Hypothesis Density，PHD)滤波器。
PHD滤波方法将集函数积分方法变换为单个目标的积分，它首先跟踪整个目标群，随后再去检测每个变量，然而PHD滤波也存在一些问题，如漏检敏感，无数目分布等。针对这一问题，Ronald P.S.Mahler教授提出了势概率假设密度滤波器(CPHD)，相比较PHD滤波器，CPHD滤波器能在传递概率密度假设函数的基础上对目标的数目分布也进行更新，在目标的估计方面做的比PHD更好。由于高斯混合势概率假设密度(GM-CPHD)滤波器假设的目标运动符合线性模型，而在实际雷达量测信息中，数据给与的三个量斜距、方位角、高地角在雷达坐标系中是非线性，这就需要对雷达量测信息进行量测转换。
虽然CPHD滤波器的性能优于PHD滤波器，但是CPHD滤波器的计算量也比PHD滤波器多得多，在CPHD滤波器中，计算复杂度为O(NM3)，而相比之下在PHD滤波器中的计算复杂度只有O(NM)，其中N为跟踪的目标数目，M为当前观测集中的观测数目。从计算复杂度可以看出，CPHD滤波器比PHD滤波器大，而当前观测集中的观测数目M作为计算复杂度的关键部分也比跟踪的目标数目N要大。
三、发明内容
本发明的目的在于提供一种基于无偏转换量测的带椭球门限高斯混合势概率假设密度滤波方法，通过无偏转换量测思想解决雷达量测信息的非线性，通过椭球门限减小计算量。
实现本发明目的的技术解决方案为：一种改进的高斯混合势概率假设密度滤波方法，包括以下步骤：
步骤1，对于多目标跟踪，目标状态集Xk＝{xk,1,…,xk,m(k)}，m(k)是目标状态向量个数，目标状态随机有限集其中Sk(Xk-1)、Nk(Xk-1)分别为原保存和新产生的目标随机有限集；k时刻新生的真实目标强度函数其中分别代表第i个新生目标的权值、均值和协方差矩阵；Jγ,k表示新生目标的总数，真实目标和杂波源的新生概率假设密度、势分布分别为
步骤2，初始化初始目标的概率假设密度D0(x)及势分布p0(n)；
步骤3，对目标状态集Xk在k+1时刻的概率假设密度及势分布进行预测，其中k≥1，得到k+1时刻的概率假设密度Dk+1|k(x)及势分布pk+1|k(n)；
步骤4，对目标状态集Xk在k+1时刻的概率假设密度及势分布进行更新，其中k≥1，得到k+1时刻的概率假设密度Dk+1(x)及势分布pk+1(n)，对真实协方差矩阵和真实偏差进行无偏转换，设置椭球门限对量测集合进行简化，减小当前观测集中的观测数目M；
步骤5，对目标的强度函数υk+1(x)的高斯项进行修剪合并，提取目标状态估计进行性能评估；
步骤6，重复步骤3～5，跟踪目标直到目标消失。
进一步地，步骤2所述初始化初始目标的概率假设密度D0(x)及势分布p0(n)，具体如下：
初始目标的概率假设密度D0(x)及势分布p0(n)的关系为：
D0(x)＝n0·s0(x)
其中，s0(x)为一概率密度，s0(x)峰值对应先验的目标位置；p0(n)是目标数n的概率分布，p0(n)期望值为n0，即：
在高斯势概率假设密度方法中，初始概率假设密度D0(x)符合高斯分布，D0(x)由每个目标的正态分布概率和表示；而初始势分布选择为二项分布，则：
其中，[0,L]为目标满足均匀分布的区间，n0为初始目标数的推测值，n0＝Lq0，q0为二项式分布发生概率，为伯努利分布、CL,n为分布系数。
进一步地，步骤3所述对目标状态集Xk在k+1时刻的概率假设密度及势分布进行预测，其中k≥1，得到k+1时刻的概率假设密度Dk+1|k(x)及势分布pk+1|k(n)，具体如下：在k时刻，已知的参数有概率假设密度Dk(x)、目标数的期望nk、势分布pk(x)，k-1时刻存活下来的目标状态集Xk的概率假设密度Dk+1|k(ξ)为：
Dk+1|k(ξ)＝∫ps(x')·fk+1|k(x|x')·Dk|k(x')dx'
其中，k≥1，ps(x')表示目标存活概率；fk+1|k(x|x')表示单目标马儿可夫转移密度；Dk|k(x')表示前一时刻目标状态集Xk的概率假设密度；
则k时刻目标状态集Xk的概率假设密度Dk+1|k(x)为：
Dk+1|k(x)＝b(x)+∫ps(x')·fk+1|k(x|x')·Dk|k(x')dx'
势分布pk+1|k(n)如下：
其中，代表第j个目标的权值，Nmax代表势分布的最大可能数，pk(l)代表前一时刻即k时刻的目标存活概率，代表二项式系数；
目标数的期望值预测如下：
其中： 分别代表期望的新生目标数和存活目标数。
进一步地，步骤4所述对目标状态集Xk在k+1时刻的概率假设密度及势分布进行更新，其中k≥1，得到k+1时刻的概率假设密度Dk+1(x)及势分布pk+1(n)，对真实协方差矩阵和真实偏差进行无偏转换，设置椭球门限对量测集合进行简化，减小当前观测集中的观测数目M，具体如下：
(4.1)在预测目标状态的强度函数和预测目标状态集的势分布pk+1|k已知且高斯混合的情况下，得到CPHD滤波器的更新方程；
势分布更新：
其中ωk+1|k为单个高斯分量的权值、Zk为k时刻目标量测值、pk+1|k为预测的势分布、为一个与权值和量测值有关的中间变量；
目标数更新：
其中，为k时刻目标数目，为k+1时刻的目标数目预测值，ωk+1|k为单个高斯分量的权值，Zk为k时刻目标量测值，pk+1|k为预测的势分布，为一个与权值和量测值有关的中间变量，pD,k+1为目标检测概率，H为量测矩阵，为目标状态协方差矩阵，为GM-CPHD函数中单个高斯分量的均值，为量测偏差，R为量测协方差；
(4.2)真实协方差矩阵和真实偏差中添加偏差补偿因子，进行无偏估计，得到改进的和代入更新方程；
μ＝E[υm|rm,θm,βm]＝[μx,μy,μz]T
其中，E为期望函数，υm为高斯混合预测强度函数，rm为目标真实斜距，θm、βm分别为高地角和方位角，μx、μy、μz分别为偏差μ在x、y、z轴上的投影，λθ、λη分别为高地角和方位角方向上的偏差补偿因子，分别为斜距、高地角和方位角的测量值，等为协方差矩阵中的量；
(4.3)设置椭球门限去掉所有和已知被探测目标不相关的量测数据，设γ为椭球跟踪门限的大小，观测维数M，残差协方差矩阵S(k)，则残差向量d(k)的范数为g(k)＝dT(k)S-1(k)d(k)，g(k)为服从自由度为M的分布，当g(k)≤γ时，椭球跟踪门规则满足。
进一步地，步骤5所述对目标的强度函数υk+1(x)的高斯项进行修剪合并，提取目标状态估计进行性能评估，具体如下：
(5.1)将小于权值τ的高斯成分滤掉，
式中为单个高斯分量的权值，为GM-CPHD函数中单个高斯分量的均值，为目标状态协方差矩阵，x为目标状态集，N为与目标状态、均值和协方差矩阵有关的函数，为k+1时刻的概率假设密度函数；
(5.2)当高斯成分间的距离小于阈值U时，将这些高斯成分合并；
(5.3)最后状态估计，目标的状态估计是提取权值大于τ1的高斯成分，其提取公式如下：其中，为单个高斯分量的权值，为GM-CPHD函数中单个高斯分量的均值，τ1为设定的门限。
本发明与现有技术相比，其显著优点在于：(1)在传递PHD函数的同时传递目标数分布，并对概率假设密度及势分布进行预测和更新，可以在杂波环境下对目标的状态和数目准确的估计；(2)对真实协方差矩阵和真实偏差进行无偏转换，可以克服GM-CPHD滤波器在雷达量测信息的非线性问题；(3)椭球门限的设定，对减小滤波器的计算量起到良好的作用，使GM-CPHD滤波器的工程应用成为可能。
四、附图说明
图1是本发明改进的高斯混合势概率假设密度滤波方法的流程图。
图2是本发明中雷达与目标相对位置示意图。
图3是本发明实施例中目标量测信息图。
图4是本发明实施例中的滤波结果图。
图5是本发明实施例与传统方法估计的目标个数图。
图6是本发明实施例中PHD和CPHD滤波器的OSPA距离对比对比图。
图7是传统CPHD滤波器和本发明算法运行耗费时间对比图。
五、具体实施方式
结合图1，本发明改进的高斯混合势概率假设密度滤波方法，包括以下步骤：
步骤1，对于多目标跟踪，目标状态集Xk＝{xk,1,…,xk,m(k)}，m(k)是目标状态向量个数，目标状态随机有限集其中Sk(Xk-1)、Nk(Xk-1)分别为原保存和新产生的目标随机有限集；k时刻新生的真实目标强度函数其中分别代表第i个新生目标的权值、均值和协方差矩阵；表示新生目标的总数，真实目标和杂波源的新生概率假设密度、势分布分别为
步骤2，初始化初始目标的概率假设密度D0(x)及势分布p0(n)，具体如下：
初始目标的概率假设密度D0(x)及势分布p0(n)的关系为：
D0(x)＝n0·s0(x)
其中，s0(x)为一概率密度，s0(x)峰值对应先验的目标位置；p0(n)是目标数n的概率分布，p0(n)期望值为n0，即：
在高斯势概率假设密度方法中，初始概率假设密度D0(x)符合高斯分布，D0(x)由每个目标的正态分布概率和表示；而初始势分布选择为二项分布，则：
其中，[0,L]为目标满足均匀分布的区间，n0为初始目标数的推测值，n0＝Lq0，q0为二项式分布发生概率，为伯努利分布、CL,n为分布系数。
步骤3，对目标状态集Xk在k+1时刻的概率假设密度及势分布进行预测，其中k≥1，得到k+1时刻的概率假设密度Dk+1|k(x)及势分布pk+1|k(n)，具体如下：
在势概率假设密度滤波器中：目标运动是独立的、不相关的，雷达与目标相对位置如图2所示，即目标x在k时刻雷达中出现的概率为bk(x)是确定的，与目标个数、状态等无关，同理目标消失的概率也一样。
在k时刻，已知的参数有：概率假设密度Dk(x)、目标数的期望nk、势分布pk(x)，k-1时刻存活下来的目标状态集Xk的概率假设密度Dk+1|k(ξ)为：
Dk+1|k(ξ)＝∫ps(x')·fk+1|k(x|x')·Dk|k(x')dx'
其中，k≥1，ps(x')表示目标存活概率，一般取0.9；fk+1|k(x|x')表示单目标马儿可夫转移密度；Dk|k(x')表示前一时刻目标状态集Xk的概率假设密度；
则k时刻目标状态集Xk的概率假设密度Dk+1|k(x)为：
Dk+1|k(x)＝b(x)+∫ps(x')·fk+1|k(x|x')·Dk|k(x')dx'
势分布pk+1|k(n)如下：
其中，代表第j个目标的权值，Nmax代表势分布的最大可能数，pk(l)代表前一时刻即k时刻的目标存活概率，代表二项式系数；
目标数的期望值预测如下：
其中： 分别代表期望的新生目标数和存活目标数。
步骤4，对目标状态集Xk在k+1时刻的概率假设密度及势分布进行更新，其中k≥1，得到k+1时刻的概率假设密度Dk+1(x)及势分布pk+1(n)，对真实协方差矩阵和真实偏差进行无偏转换，设置椭球门限对量测集合进行简化，减小当前观测集中的观测数目M，具体如下：
(4.1)在预测目标状态的强度函数和预测目标状态集的势分布pk+1|k已知且高斯混合的情况下，得到CPHD滤波器的更新方程；
势分布更新：
其中ωk+1|k为单个高斯分量的权值、Zk为k时刻目标量测值、pk+1|k为预测的势分布、为一个与权值和量测值有关的中间变量；
目标状态的强度函数更新：
目标数更新：
其中
其中，δj(·)为均衡函数，κk(·)为杂波强度函数，为k时刻目标数目，为k+1时刻的目标数目预测值，ωk+1|k为单个高斯分量的权值，Zk为k时刻目标量测值，pk+1|k为预测的势分布，为一个与权值和量测值有关的中间变量，pD,k+1为目标检测概率，H为量测矩阵，为目标状态协方差矩阵，为GM-CPHD函数中单个高斯分量的均值，为量测偏差，R为量测协方差。
(4.2)由于GM-CPHD假设目标的运动模型服从线性高斯，协方差矩阵Rk+1内的元素是根据传感器特性设定好的定值。而接收到的目标信息都是在瞄具极坐标系下的量测信息。瞄具下的目标信息为径向距离、方位角和高低角，而这三个信息在瞄具坐标系下是非线性的，为了将GM-CPHD算法运用到瞄具坐标系下的多目标跟踪，将无偏转换思想融入到传统的GM-CPHD算法中，以求对目标状态估计效果达到最佳。
真实协方差矩阵和真实偏差中添加偏差补偿因子，进行无偏估计，得到改进的和代入更新方程：
μ＝E[υm|rm,θm,βm]＝[μx,μy,μz]T
其中，
分别是和的协方差。
其中，E为期望函数，υm为高斯混合预测强度函数，rm为目标真实斜距，θm、βm分别为高地角和方位角，μx、μy、μz分别为偏差μ在x、y、z轴上的投影，λθ、λη分别为高地角和方位角方向上的偏差补偿因子，分别为斜距、高地角和方位角的测量值，等为协方差矩阵中的量。
(4.3)在CPHD滤波器中，算法的计算量主要来自每个更新周期都要计算M+1个元素的均衡函数(即前文提到的δj(·))，其计算复杂度为O(NM3)，而相比之下在PHD滤波器中的计算复杂度只有O(NM)，其中N为跟踪的目标数目，M为当前观测集中的观测数目，由复杂度公式可知，减小M值能更有效的降低计算复杂度，而采用椭球门限不失为一种不错的方法。
设置椭球门限去掉所有和已知被探测目标不相关的量测数据，设γ为椭球跟踪门限的大小，观测维数M，残差协方差矩阵S(k)，则残差向量d(k)的范数为g(k)＝dT(k)S-1(k)d(k)，g(k)为服从自由度为M的分布，当g(k)≤γ时，椭球跟踪门规则满足。
首先计算出目标量测数据落入椭球跟踪门内的概率PG：
结合椭球跟踪门规则可得：
其中：
椭球跟踪门的体积为：
其中，系数由观测空间的维数nz决定(c1＝2，c2＝π，c3＝4π/3)。
然后设置椭球门限减小M值，本发明中的观测空间的维数nz＝3，此时门限γ与目标量测数据落入椭球跟踪门内的概率PG有如下关系：
一般取pG＝0.99，计算得到此时γ＝11.34。
椭球门限所包含的区域则落入联合门限区域的量测集合就可以表示为：与传统算法相比，采用椭球门限可能会去掉所有和已知被探测目标不相关的量测数据，减小计算量。
最后更新目标状态的强度函数υk+1|k(x)及势分布pk+1|k(n)，由于量测集合由Zk变化为引起杂波强度函数变化为则更新公式中
杂波强度的变化直接影响到势概率假设密度滤波器中高斯分量对应权值和量测状态集的变化，而后者的大小M恰恰是CPHD滤波器计算复杂度O(NM3)的关键部分，所以通过椭球门限可以有效的减小滤波器的计算量。
步骤5，对目标的强度函数υk+1(x)的高斯项进行修剪合并，提取目标状态估计进行性能评估。
(5.1)首先修剪，将小于权值τ的高斯成分滤掉，
式中为单个高斯分量的权值，为GM-CPHD函数中单个高斯分量的均值，为目标状态协方差矩阵，x为目标状态集，N为与目标状态、均值和协方差矩阵有关的函数，为k+1时刻的概率假设密度函数；
(5.2)然后合并，当高斯成分间的距离小于阈值U时，将这些高斯成分合并。设置l＝0,且重复以下步骤：
l＝l+1
I＝I\L，直到I＝φ。
如果l≤Jmax，求取的便是合并后的高斯成分。否则需要用Jmax个高斯成分代替这些代替后的高斯成分的权值是合并后最大的前Jmax个。
(5.3)最后状态估计，目标的状态估计是提取权值大于τ1的高斯成分，其提取公式如下：其中，为单个高斯分量的权值，为GM-CPHD函数中单个高斯分量的均值，τ1为设定的门限。采用最大后验(MAP)估计器可以有效的忽略虚警等，锁定稳定和精确的主模式，
性能评价：对于多目标跟踪算法的性能评价指标，一般采用均方误差(MSE)、均方根误差(RMSE)、圆丢失概率(CPEP)、Wasserstein距离和OSPA距离等。一般OSPA距离可以通过水平参数c来很好地体现多目标跟踪算法对位置和目标数目跟踪的性能。OSPA距离定义如下：
其中，
步骤6，重复步骤3～5，跟踪目标直到目标消失。
实施例1
本发明的效果可以通过以下仿真实验进一步说明：
1、仿真条件
假设目标状态其中位置单位是m，速度单位是m/s。本仿真有四个目标，各个目标运动模型是CA模型，且四个目标的初始状态如下：
Xt1＝[0m,0m,0m,-3m/s,-3m/s,-3m/s,-0.5m/s2,0m/s2,0m/s2]T；
Xt2＝[0m,0m,0m,-3m/s,3m/s,3m/s,-0.3m/s2,0m/s2,0m/s2]T；
Xt3＝[0m,0m,0m,3m/s,-3m/s,-3m/s,0.5m/s2,0m/s2,0m/s2]T；
Xt4＝[0m,0m,0m,-3m/s,-3m/s,3m/s,-0.8m/s2,0m/s2,0m/s2]T；
目标的运动方程Xk＝FXk-1+wk，其中存活目标转移矩阵F：
仿真实验假设目标1的出生时刻为第0s，死亡时刻为第7s，目标2的出生时刻为第7s，死亡时刻为25s，目标3的出生时刻为第11s，死亡时刻为37s，目标4的出生时刻为25s，死亡时刻为第40s。设雷达采样周期T＝1s，径向距离量测方差方位角与高低角量测方差取检测概率PD＝0.99，目标存活概率PS＝0.9，合并阈值U＝4，裁剪阈值τ＝1e-5，状态估计阈值τ1＝0.5，最大高斯数Jmax＝100。
新出现的目标符合泊松过程：
bk(x)＝N(x,mb,Pb)
mb＝[0m,0m,0m,0m/s,0m/s,0m/s,0m/s2,0m/s2,0m/s2]T
Pb＝diag(104×[5,5,5,1,1,1,0.1,0.1,0.1])，Q＝diag(10-2×[1,1,1])
2、仿真内容和结果分析
生成的目标量测信息如图3所示，量测信息中含有目标信息和杂波信息。
目标估计值与目标真实值对比如图4所示，可以看出目标估计点在真实值附近波动。其中黑色为基于无偏转换的带椭球门限高斯混合概率假设密度滤波器估计，红色为本发明(基于无偏转换的带椭球门限高斯混合势概率假设密度滤波器)估计值，从图中可以看出，对于目标状态估计，本发明比前者更准确。
算法各个时刻估计的目标个数结果如图5所示，可以看出，该算法对于目标数目的估计非常准确，比PHD滤波器做的更好。
图6为PHD和CPHD滤波器的OSPA距离对比，在多目标跟踪算法中，OSPA距离是衡量算法性能指标的重要参数。从图中看出，本发明性能优于PHD滤波器。通过多次实验得到传统CPHD滤波器和本发明算法运行耗费时间对比，由图7可知，相比较于传统方法，本发明算法运行时间减小5％左右，具有更好的工程应用前景。
Claims (5)
1.一种改进的高斯混合势概率假设密度滤波方法，其特征在于，包括以下步骤：
步骤1，对于多目标跟踪，目标状态集Xk＝{xk,1,…,xk,m(k)}，m(k)是目标状态向量个数，目标状态随机有限集Ξk＝Sk(Xk-1)∪Nk(Xk-1)，其中Sk(Xk-1)、Nk(Xk-1)分别为原保存和新产生的目标随机有限集；k时刻新生的真实目标强度函数其中分别代表第i个新生目标的权值、均值和协方差矩阵；Jγ,k表示新生目标的总数，真实目标和杂波源的新生概率假设密度、势分布分别为
步骤2，初始化初始目标的概率假设密度D0(x)及势分布p0(n)；
步骤3，对目标状态集Xk在k+1时刻的概率假设密度及势分布进行预测，其中k≥1，得到k+1时刻的概率假设密度Dk+1|k(x)及势分布pk+1|k(n)；
步骤4，对目标状态集Xk在k+1时刻的概率假设密度及势分布进行更新，其中k≥1，得到k+1时刻的概率假设密度Dk+1(x)及势分布pk+1(n)，对真实协方差矩阵和真实偏差进行无偏转换，设置椭球门限对量测集合进行简化，减小当前观测集中的观测数目M；
步骤5，对目标的强度函数υk+1(x)的高斯项进行修剪合并，提取目标状态估计进行性能评估；
步骤6，重复步骤3～5，跟踪目标直到目标消失。
2.根据权利要求1所述的改进的高斯混合势概率假设密度滤波方法，其特征在于，步骤2所述初始化初始目标的概率假设密度D0(x)及势分布p0(n)，具体如下：
初始目标的概率假设密度D0(x)及势分布p0(n)的关系为：
D0(x)＝n0·s0(x)
其中，s0(x)为一概率密度，s0(x)峰值对应先验的目标位置；p0(n)是目标数n的概率分布，p0(n)期望值为n0，即：
在高斯势概率假设密度方法中，初始概率假设密度D0(x)符合高斯分布，D0(x)由每个目标的正态分布概率和表示；而初始势分布选择为二项分布，则：
其中，[0,L]为目标满足均匀分布的区间，n0为初始目标数的推测值，n0＝Lq0，q0为二项式分布发生概率，为伯努利分布、CL,n为分布系数。
3.根据权利要求1所述的改进的高斯混合势概率假设密度滤波方法，其特征在于，步骤3所述对目标状态集Xk在k+1时刻的概率假设密度及势分布进行预测，其中k≥1，得到k+1时刻的概率假设密度Dk+1|k(x)及势分布pk+1|k(n)，具体如下：在k时刻，已知的参数有概率假设密度Dk(x)、目标数的期望nk、势分布pk(x)，k-1时刻存活下来的目标状态集Xk的概率假设密度Dk+1|k(ξ)为：
Dk+1|k(ξ)＝∫ps(x')·fk+1|k(x|x')·Dk|k(x')dx'
其中，k≥1，ps(x')表示目标存活概率；fk+1|k(x|x')表示单目标马儿可夫转移密度；Dk|k(x')表示前一时刻目标状态集Xk的概率假设密度；
则k时刻目标状态集Xk的概率假设密度Dk+1|k(x)为：
Dk+1|k(x)＝b(x)+∫ps(x')·fk+1|k(x|x')·Dk|k(x')dx'
势分布pk+1|k(n)如下：
其中，代表第j个目标的权值，Nmax代表势分布的最大可能数，pk(l)代表前一时刻即k时刻的目标存活概率，代表二项式系数；
目标数的期望值预测如下：
其中： 分别代表期望的新生目标数和存活目标数。
4.根据权利要求1所述的改进的高斯混合势概率假设密度滤波方法，其特征在于，步骤4所述对目标状态集Xk在k+1时刻的概率假设密度及势分布进行更新，其中k≥1，得到k+1时刻的概率假设密度Dk+1(x)及势分布pk+1(n)，对真实协方差矩阵和真实偏差进行无偏转换，设置椭球门限对量测集合进行简化，减小当前观测集中的观测数目M，具体如下：
(4.1)在预测目标状态的强度函数和预测目标状态集的势分布pk+1|k已知且高斯混合的情况下，得到CPHD滤波器的更新方程；
势分布更新：
其中ωk+1|k为单个高斯分量的权值、Zk为k时刻目标量测值、pk+1|k为预测的势分布、为一个与权值和量测值有关的中间变量；
目标数更新：
其中，为k时刻目标数目，为k+1时刻的目标数目预测值，ωk+1|k为单个高斯分量的权值，Zk为k时刻目标量测值，pk+1|k为预测的势分布，为一个与权值和量测值有关的中间变量，pD,k+1为目标检测概率，H为量测矩阵，为目标状态协方差矩阵，为GM-CPHD函数中单个高斯分量的均值，为量测偏差，R为量测协方差；
(4.2)真实协方差矩阵和真实偏差中添加偏差补偿因子，进行无偏估计，得到改进的和代入更新方程；
μ＝E[υm|rm,θm,βm]＝[μx,μy,μz]T
其中，E为期望函数，υm为高斯混合预测强度函数，rm为目标真实斜距，θm、βm分别为高地角和方位角，μx、μy、μz分别为偏差μ在x、y、z轴上的投影，λθ、λη分别为高地角和方位角方向上的偏差补偿因子，分别为斜距、高地角和方位角的测量值，等为协方差矩阵中的量；
(4.3)设置椭球门限去掉所有和已知被探测目标不相关的量测数据，设γ为椭球跟踪门限的大小，观测维数M，残差协方差矩阵S(k)，则残差向量d(k)的范数为g(k)＝dT(k)S-1(k)d(k)，g(k)为服从自由度为M的分布，当g(k)≤γ时，椭球跟踪门规则满足。
5.根据权利要求1所述的改进的高斯混合势概率假设密度滤波方法，其特征在于，步骤5所述对目标的强度函数υk+1(x)的高斯项进行修剪合并，提取目标状态估计进行性能评估，具体如下：
(5.1)将小于权值τ的高斯成分滤掉，
式中为单个高斯分量的权值，为GM-CPHD函数中单个高斯分量的均值，为目标状态协方差矩阵，x为目标状态集，N为与目标状态、均值和协方差矩阵有关的函数，为k+1时刻的概率假设密度函数；
(5.2)当高斯成分间的距离小于阈值U时，将这些高斯成分合并；
(5.3)最后状态估计，目标的状态估计是提取权值大于τ1的高斯成分，其提取公式如下：其中，为单个高斯分量的权值，为GM-CPHD函数中单个高斯分量的均值，τ1为设定的门限。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610325272.8A CN106022340A (zh) | 2016-05-17 | 2016-05-17 | 一种改进的高斯混合势概率假设密度滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610325272.8A CN106022340A (zh) | 2016-05-17 | 2016-05-17 | 一种改进的高斯混合势概率假设密度滤波方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN106022340A true CN106022340A (zh) | 2016-10-12 |
Family
ID=57097208
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610325272.8A Pending CN106022340A (zh) | 2016-05-17 | 2016-05-17 | 一种改进的高斯混合势概率假设密度滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106022340A (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106802414A (zh) * | 2016-12-19 | 2017-06-06 | 姜秋喜 | 基于高斯滤波的机动目标跟踪方法 |
| CN107831490A (zh) * | 2017-12-01 | 2018-03-23 | 南京理工大学 | 一种改进的多扩展目标跟踪方法 |
| CN108227750A (zh) * | 2017-12-20 | 2018-06-29 | 西安石油大学 | 一种地面目标实时跟踪性能评估方法及系统 |
| CN108717702A (zh) * | 2018-04-24 | 2018-10-30 | 西安工程大学 | 基于分段rts的概率假设密度滤波平滑方法 |
| CN108732564A (zh) * | 2018-03-13 | 2018-11-02 | 中国电子科技集团公司第二十八研究所 | 一种双雷达修正序贯高斯混合概率假设密度滤波方法 |
| CN109034356A (zh) * | 2018-07-23 | 2018-12-18 | 北京理工大学 | 基于最近邻域法关联和高斯波束体积的昆虫密度统计方法 |
| CN109946694A (zh) * | 2019-03-22 | 2019-06-28 | 哈尔滨工业大学 | 基于随机有限集的圆周sar多目标跟踪方法 |
| CN110334472A (zh) * | 2019-07-15 | 2019-10-15 | 中国人民解放军国防科技大学 | 一种群运动趋势辅助的集势概率假设密度滤波方法 |
| CN111736145A (zh) * | 2020-06-28 | 2020-10-02 | 电子科技大学 | 一种基于高斯混合概率假设密度滤波的多机动目标多普勒雷达跟踪方法 |
| CN112654979A (zh) * | 2020-04-29 | 2021-04-13 | 华为技术有限公司 | 数据关联方法与装置 |
| CN113156279A (zh) * | 2021-04-09 | 2021-07-23 | 江苏大学 | 一种基于概率假设密度滤波器的配电柜局部放电定位方法 |
-
2016
- 2016-05-17 CN CN201610325272.8A patent/CN106022340A/zh active Pending
Cited By (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106802414B (zh) * | 2016-12-19 | 2019-07-12 | 姜秋喜 | 基于高斯滤波的机动目标跟踪方法 |
| CN106802414A (zh) * | 2016-12-19 | 2017-06-06 | 姜秋喜 | 基于高斯滤波的机动目标跟踪方法 |
| CN107831490A (zh) * | 2017-12-01 | 2018-03-23 | 南京理工大学 | 一种改进的多扩展目标跟踪方法 |
| CN108227750B (zh) * | 2017-12-20 | 2021-02-05 | 西安石油大学 | 一种地面目标实时跟踪性能评估方法及系统 |
| CN108227750A (zh) * | 2017-12-20 | 2018-06-29 | 西安石油大学 | 一种地面目标实时跟踪性能评估方法及系统 |
| CN108732564A (zh) * | 2018-03-13 | 2018-11-02 | 中国电子科技集团公司第二十八研究所 | 一种双雷达修正序贯高斯混合概率假设密度滤波方法 |
| CN108732564B (zh) * | 2018-03-13 | 2020-04-17 | 中国电子科技集团公司第二十八研究所 | 一种双雷达修正序贯高斯混合概率假设密度滤波方法 |
| CN108717702A (zh) * | 2018-04-24 | 2018-10-30 | 西安工程大学 | 基于分段rts的概率假设密度滤波平滑方法 |
| CN109034356A (zh) * | 2018-07-23 | 2018-12-18 | 北京理工大学 | 基于最近邻域法关联和高斯波束体积的昆虫密度统计方法 |
| CN109034356B (zh) * | 2018-07-23 | 2021-06-01 | 北京理工大学 | 基于最近邻域法关联和高斯波束体积的昆虫密度统计方法 |
| CN109946694A (zh) * | 2019-03-22 | 2019-06-28 | 哈尔滨工业大学 | 基于随机有限集的圆周sar多目标跟踪方法 |
| CN110334472A (zh) * | 2019-07-15 | 2019-10-15 | 中国人民解放军国防科技大学 | 一种群运动趋势辅助的集势概率假设密度滤波方法 |
| CN112654979A (zh) * | 2020-04-29 | 2021-04-13 | 华为技术有限公司 | 数据关联方法与装置 |
| CN111736145A (zh) * | 2020-06-28 | 2020-10-02 | 电子科技大学 | 一种基于高斯混合概率假设密度滤波的多机动目标多普勒雷达跟踪方法 |
| CN111736145B (zh) * | 2020-06-28 | 2022-04-19 | 电子科技大学 | 一种基于高斯混合概率假设密度滤波的多机动目标多普勒雷达跟踪方法 |
| CN113156279A (zh) * | 2021-04-09 | 2021-07-23 | 江苏大学 | 一种基于概率假设密度滤波器的配电柜局部放电定位方法 |
| CN113156279B (zh) * | 2021-04-09 | 2023-01-17 | 江苏大学 | 一种基于概率假设密度滤波器的配电柜局部放电定位方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106022340A (zh) | 一种改进的高斯混合势概率假设密度滤波方法 | |
| CN106407677B (zh) | 一种测量数据丢失情况下的多目标跟踪方法 | |
| Li et al. | Kalman filter and its application | |
| CN110503071B (zh) | 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法 | |
| CN106487358B (zh) | 一种机动目标转弯跟踪方法 | |
| US20230297642A1 (en) | Bearings-only target tracking method based on pseudo-linear maximum correlation entropy kalman filtering | |
| CN107831490A (zh) | 一种改进的多扩展目标跟踪方法 | |
| CN108896986B (zh) | 一种基于预测值的量测转换序贯滤波机动目标跟踪方法 | |
| CN106372646B (zh) | 基于srck-gmcphd滤波的多目标跟踪方法 | |
| CN105405151A (zh) | 基于粒子滤波和加权Surf的抗遮挡目标跟踪方法 | |
| CN101975575A (zh) | 基于粒子滤波的被动传感器多目标跟踪方法 | |
| CN111722214A (zh) | 雷达多目标跟踪phd实现方法 | |
| CN110889862A (zh) | 一种网络传输攻击环境中多目标跟踪的组合测量方法 | |
| CN115204212A (zh) | 一种基于stm-pmbm滤波算法的多目标跟踪方法 | |
| Petetin et al. | Marginalized particle PHD filters for multiple object Bayesian filtering | |
| CN112379350A (zh) | 智能车辆毫米波雷达多目标跟踪方法、装置及设备 | |
| CN109375160B (zh) | 纯方位无源定位中一种测角误差估计方法 | |
| Williams | Experiments with graphical model implementations of multiple target multiple Bernoulli filters | |
| García-Fernández et al. | A time-weighted metric for sets of trajectories to assess multi-object tracking algorithms | |
| CN114445456B (zh) | 基于部分模型的数据驱动智能机动目标跟踪方法及装置 | |
| CN115544425A (zh) | 一种基于目标信噪比特征估计的鲁棒多目标跟踪方法 | |
| CN115619825A (zh) | 地面多目标跟踪状态及轨迹确定方法 | |
| CN114611068A (zh) | 一种高机动目标跟踪方法 | |
| Yan et al. | Multi-level Deep Learning Kalman Filter | |
| Wang et al. | State Estimation under Outliers by the Maximum Correntropy Extended Kalman Filter |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20161012 |