CN115619825A - 地面多目标跟踪状态及轨迹确定方法 - Google Patents

地面多目标跟踪状态及轨迹确定方法 Download PDF

Info

Publication number
CN115619825A
CN115619825A CN202211227444.XA CN202211227444A CN115619825A CN 115619825 A CN115619825 A CN 115619825A CN 202211227444 A CN202211227444 A CN 202211227444A CN 115619825 A CN115619825 A CN 115619825A
Authority
CN
China
Prior art keywords
target
gaussian
measurement
label
jth
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
Application number
CN202211227444.XA
Other languages
English (en)
Inventor
母晓慧
罗飞
崔雷
郝雅雯
谢红宇
张晓璐
张超蕊
董琦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kirin Software Co Ltd
Original Assignee
Kirin Software Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Kirin Software Co Ltd filed Critical Kirin Software Co Ltd
Priority to CN202211227444.XA priority Critical patent/CN115619825A/zh
Publication of CN115619825A publication Critical patent/CN115619825A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/246Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/277Analysis of motion involving stochastic approaches, e.g. using Kalman filters
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20024Filtering details

Landscapes

  • Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

地面多目标跟踪状态及轨迹确定方法,包括如下步骤:采用GM来实现δ‑GLMB,将用于传递的伯努利参数利用多个高斯项加权求和的方式来代替;通过非线性滤波器CKF对GM‑δ‑GLMB中预测空间分布密度函数的存活部分所分解的高斯项进行预测;获取上一时刻的量测与当前时刻量测集中的统计距离,提取与预测量测相关的量测集,依据统计距离的大小对阈值内的量测赋予不同的权重,合并不同权重的量测;利用CKF对δ‑GLMB的高斯分量进行更新;对高斯项进行剪枝和合并;对多个目标状态进行提取并估计目标运动轨迹。本发明在强杂波密度场景下,能够实现多目标实时跟踪,在计算效率方面有明显的提升。

Description

地面多目标跟踪状态及轨迹确定方法
技术领域
本发明涉及多目标轨迹预测方法,具体涉及地面多目标跟踪状态及轨迹确定方法。
背景技术
1994年,美国学者Ronald P.S.Mahler在RFS框架下提出了基于随机集理论的多目标跟踪方法,该类算法在随机有限集统一框架下,将目标状态与量测值建模为随机有限集的形式,将多目标密度建模为集合密度,避免了经典多目标跟踪方法显式的数据关联,较好避免组合爆炸问题。
在2003年,Ba-Tuong Vo等人提出了概率假设密度(Probability HypothesisDensity,PHD)滤波,该滤波迭代传递多目标的一阶矩近似多目标后验密度,传播目标的强度函数,设定目标的数目服从泊松分布,对目标的运动状态及数目实时估计。势概率假设密度(Cardinality Probability Hypothesis Density,CPHD)是依据PHD算法加入了目标数目分布的传递,解决了PHD在目标跟踪中对传感器量测漏检敏感的问题,可以提高目标跟踪滤波算法的精度。与PHD与CPHD算法不同的是多目标多伯努利(Multi-Target Multi-Bernoulli Filter,MeMBer)算法是递归传递近似多目标后验密度的多伯努利分布随机有限集参数(存在概率与概率密度函数)进行多目标后验密度的实时估计,该算法不需要进行聚类运算,能够有效减低计算复杂度因而备受广大学者关注,为基于RFS理论的多目标跟踪方法提供了一种新思路。MeMBer算法的跟踪性能优于PHD,而计算复杂度比CPHD小。
针对MeMBer算法中进行势估计时存在目标数目过估的问题,Vo等首次提出了势平衡的多目标多伯努利滤波(Cardinality Balanced Multi-target Multi-Bernoulli,CBMeMBer)算法,此算法求得概率密度函数的近似解是通过平衡目标的势偏差。但是MeMBer与CBMeMBer算法仅能估计目标状态与数目,不能区分目标跟踪轨迹。为此,Ba-Tuong Vo等人将MHT思想和随机有限集(RFS) 理论相结合,提出δ-广义标签多伯努利(δ-GLMB)滤波算法。该滤波算法在跟踪精度方面比其他基于随机有限集的滤波算法要优,而且能够自发生成目标轨迹,是第一个具有理论证明的有精确闭式解的多目标跟踪滤波器。
由于δ-GLMB在滤波过程中存在集合积分运算,很难求得解析解,可通过GM与非线性滤波器(如卡尔曼滤波、扩展卡尔曼滤波、无迹卡尔曼滤波等)结合,可以实现非线性场景下的地面目标跟踪,并可提高目标跟踪精度。在城市环境中受建筑物发射、散射等因素的影响,杂波密集,具有极大的跟踪算法的运算复杂度。为此,本发明在保障跟踪要求的同时,通过牺牲部分精度,采用容积卡尔曼(CKF)对改进的GM-δ-GLMB滤波算法的高斯参量进行预测与更新,以适应强杂波密度场景下的地面目标跟踪。
在现有的技术中,一种强杂波密度环境下的多目标跟踪算法与本发明专利同属于一个领域,该算法在概率假设密度平滑器基础上,在预测步通过量测与预测值之间的残差值确定椭圆门限值,从而获得与目标真实状态相近的量测值。
最接近的现有技术中,虽然能够实现多目标跟踪,但是该技术通过概率假设密度滤波算法做多目标跟踪,在预测和更新的具体实现过程中,目标状态提取时需要进行粒子数目聚类,而聚类算法具有不稳定和计算量大的缺陷。
发明内容
为解决已有技术存在的不足,本发明提供了一种地面多目标跟踪状态及轨迹确定方法,包括如下步骤:
步骤S1:K=0时刻,初始化参数,获取地面多机动目标状态均值、多目标状态协方差矩阵、多目标非线性运动模型及杂波参数,得到初始分量集;采用GM来实现δ-GLMB,将用于传递的伯努利参数利用多个高斯项加权求和的方式来代替;
步骤S2:K大于或等于1时,通过非线性滤波器CKF对GM-δ-GLMB中预测空间分布密度函数的存活部分所分解的高斯项进行预测,完成对初始化的地面多目标状态集合信息的预测;
步骤S3:获取步骤S2中上一时刻的量测与当前时刻量测集中的统计距离,以此为依据通过查看χ2分布表设置阈值来提取与预测量测相关的量测集,然后,依据统计距离的大小对阈值内的量测赋予不同的权重,最后合并不同权重的量测,从而基于距离加权的量测合并策略获取观测的地面目标更新步骤对应的多个传感器量测信息;
步骤S4:结合步骤S3所获取的传感器量测信息,利用CKF对δ-GLMB 的高斯分量进行更新,得到本时刻最终输出的状态集合,对多个目标状态进行提取并估计目标运动轨迹,完成一次迭代;输出下一时刻多个传感器的量测合集,转至步骤S2进行下一次迭代,直至无关于下一刻的输入;
步骤S5:对高斯项进行剪枝和合并;
步骤S6:获取各个时刻的目标状态及目标运动轨迹,对多个目标状态进行提取并估计目标运动轨迹。
其中,所述步骤S1中,采用GM来实现δ-GLMB时,先验概率密度中标签l对应的伯努利参数的GM形式为:
Figure RE-GDA0003988893070000041
其中,k为观测时间,
Figure RE-GDA0003988893070000042
表示先验概率密度中标签l对应参数所分解成的高斯项个数,j表示先验概率密度中标签l对应参数所分解成当前的高斯项,
Figure RE-GDA0003988893070000043
表示k-1时刻的第j个目标的高斯项状态均值,
Figure RE-GDA0003988893070000044
表示其对应协方差,
Figure RE-GDA0003988893070000045
表示k-1时刻标签l第j个高斯项的权重,在实际仿真中,此值为随机非负数,且满足
Figure RE-GDA0003988893070000046
Figure RE-GDA0003988893070000047
表示k-1时刻的第j个高斯项的状态均值为
Figure RE-GDA0003988893070000048
协方差为
Figure RE-GDA0003988893070000049
的高斯概率密度;y表示自变量参数。
其中,所述步骤S2包括如下步骤:
步骤S21:生成容积点集:
Figure RE-GDA0003988893070000051
其中,s=2n,n为系统状态维度,[1]=[In×n,In×n],I是n维单位矩阵,[1]t是单位矩阵的第t列向量。CKF的容积点集{ξt,wt}},利用三阶球面径向规则计算标准高斯加权积分,即
Figure RE-GDA0003988893070000052
ξt与wt都是容积点集的参数,分别为基本容积点及权值,下标t代表单位矩阵的第t列向量;
步骤S22:预测步骤,实现高斯参量
Figure RE-GDA0003988893070000053
Figure RE-GDA0003988893070000054
的传递,包括:
步骤S221,Cholesky分解:
Figure RE-GDA0003988893070000055
Figure RE-GDA0003988893070000056
是对公式(1)组成的目标轨迹中标签l在k-1时刻的第j个高斯项的协方差矩阵进行Cholesky分解后的值;
步骤S222,产生容积点:
Figure RE-GDA0003988893070000057
Figure RE-GDA0003988893070000058
表示标签l在k-1时刻的第j个高斯项的容积点;
步骤S223,计算容积点经非线性状态方程传播后的值:
Figure RE-GDA0003988893070000059
步骤S224,估计状态预测值:
Figure RE-GDA0003988893070000061
步骤S225,计算预测误差协方差矩阵:
Figure RE-GDA0003988893070000062
其中,chol(·)表示对协方差矩阵进行cholesky分解,f表示状态转移矩阵,
Figure RE-GDA0003988893070000063
表示标签l在k-1时刻的第j个高斯项的过程噪声协方差,过程噪声是相互独立的,服从均值为0、协方差矩阵为R的高斯分布。
其中,所述步骤S3中,通过如下步骤推导出单目标的量测集:
步骤S31:定义
Figure RE-GDA0003988893070000064
Figure RE-GDA0003988893070000065
的统计距离为:
Figure RE-GDA0003988893070000066
其中,
Figure RE-GDA0003988893070000067
为标签l在k-1时刻的第j个高斯项的新息协方差矩阵;
步骤S32:通过选定特定的关联阈值,利用下式的统计门提取第j个假设轨迹关联的量测:
Figure RE-GDA0003988893070000068
其中,Th为第k时刻实现预定的阈值;
步骤S33:
Figure RE-GDA0003988893070000069
其对应的权值被描述为:
Figure RE-GDA00039888930700000610
其中,
Figure RE-GDA00039888930700000611
表示单目标在k时刻标签l第j个假设轨迹的相关量测,
Figure RE-GDA00039888930700000612
表示多目标在k时刻标签l第j个假设轨迹相关量测;
Figure RE-GDA0003988893070000071
为归一化因子;
步骤S34:将量测集表示为:
Figure RE-GDA0003988893070000072
其中,
Figure RE-GDA0003988893070000073
表示k-1时刻单目标标签l第j个高斯项的传感器量测,
Figure RE-GDA0003988893070000074
表示k时刻单目标标签l第j个高斯项的传感器量测。
其中,所述步骤S31中,式(8)中的新息协方差矩阵
Figure RE-GDA0003988893070000075
通过如下式子得到:
步骤S311,更新步Cholesky分解:
Figure RE-GDA0003988893070000076
其中,
Figure RE-GDA0003988893070000077
表示对预测结果的标签l第j个高斯项的协方差矩阵进行 Cholesky分解值;
步骤S312,计算容积点:
Figure RE-GDA0003988893070000078
其中,
Figure RE-GDA0003988893070000079
为标签l在k时刻的第j个高斯项的基本容积点;
步骤S313,计算容积点经传感器量测非线性函数传播后的值:
Figure RE-GDA00039888930700000710
其中,
Figure RE-GDA0003988893070000081
指标签l在k时刻的第j个高斯项的容积点经传感器量测非线性函数传播后的值;
步骤S314,计算新息协方差:
Figure RE-GDA0003988893070000082
其中,
Figure RE-GDA0003988893070000083
为标签l在k时刻的第j个高斯项的量测噪声协方差。
其中,所述步骤S4包括:
步骤S41:更新步Cholesky分解,操作同于步骤S311;
步骤S42:计算容积点,操作同于步骤S312;
步骤S43:计算容积点经传感器量测非线性函数传播后的值,操作同于步骤S313;
步骤S44:计算信息协方差,操作同于步骤S314;
步骤S45,计算量测预测值:
Figure RE-GDA0003988893070000084
步骤S46,计算互协方差:
Figure RE-GDA0003988893070000085
步骤S47,计算卡尔曼增益:
Figure RE-GDA0003988893070000086
步骤S48,更新状态均值:
Figure RE-GDA0003988893070000091
步骤S49,更新协方差矩阵:
Figure RE-GDA0003988893070000092
至此,实现高斯参量预测步到更新步
Figure RE-GDA0003988893070000093
的转变;
其中,h(·)表示量测转移函数,
Figure RE-GDA0003988893070000094
是量测噪声协方差,量测噪声是相互独立的,服从均值为0、协方差矩阵为Q的高斯分布;
输出下一时刻的多目标量测集合,转到步骤S2进行迭代,直到无关于下一时刻的输入时,转至步骤S5。
本发明的地面多目标跟踪状态及轨迹确定方法,提出了一种改进 CKF-GM-δ-GLMB多目标跟踪滤波算法。在强杂波密度场景下,能够实现多目标实时跟踪,在计算效率方面有明显的提升,且从其运行时间随杂波密度增加的涨幅来看,本发明的运行时间涨幅略小,进一步验证了本发明的有效性。
附图说明
图1:本发明的总体思路结构图。
图2:本发明的一实施例中八个地面目标的真实轨迹图。
图3:本发明的一实施例的滤波结果图。
图4:本发明的一实施例的在x轴和y轴方向的目标滤波结果图。
图5:本发明的一实施例的OSPA距离对比图。
图6:本发明的一实施例的势估计误差对比图。
具体实施方式
为了对本发明的技术方案及有益效果有更进一步的了解,下面结合附图详细说明本发明的技术方案及其产生的有益效果。
本发明的目的是提供一种基于随机有限集的CKF-GM-δ-GLMB的改进算法,在确保多目标滤波算法滤波效果也就是满足基本精度保障的基础上,降低目标跟踪算法的计算复杂度,解决面向地面目标跟踪强杂波密度场景下的目标实时准确跟踪这一关键问题。
本发明提供一种基于随机有限集的面向强杂波密度场景的地面多目标跟踪方法,针对繁华的闹市区杂波密度高而目标个数相对较少的场景,在确保多目标滤波算法滤波效果满足基本精度保障的基础上,降低目标跟踪算法的计算复杂度,从而解决面向地面目标跟踪强杂波密度场景下的目标实时准确跟踪这一关键问题,具有良好的工程应用价值。该方法首先将传感器的初始化目标运动状态及目标数量输入到CKF-GM-δ-GLMB滤波算法的预测步骤中,基于基于距离加权的量测合并策略求得观测目标更新步骤对应的量测,然后采用CKF-GM-δ-GLMB滤波算法进行更新地面多目标的伯努利参量,从而减少多目标跟踪算法的计算复杂度,提高多目标跟踪的有效性。方法总体结构图如图1所示。
一、整体思路:本方法分为六个步骤:
(1)本方法采用GM来实现δ-GLMB,GM是将用于传递的伯努利参数利用多个高斯项加权求和的方式来代替,以解决预测与更新步骤中带有有限集积分运算而无法求得解析解的问题;K=0时刻,初始化参数,获取地面多机动目标状态均值、多目标状态协方差矩阵、多目标非线性运动模型及杂波参数,得到初始分量集;采用GM来实现δ-GLMB,将用于传递的伯努利参数利用多个高斯项加权求和的方式来代替;
(2)为了提升GM-δ-GLMB在非线性系统中的滤波效果,本方法引入非线性滤波器CKF对GM-δ-GLMB预测。GM-δ-GLMB中预测空间分布密度函数由新生部分和存活部分组成,其中存活部分的空间分布密度函数可分解成高斯项组成,然后采用CKF对δ-GLMB滤波算法先验信息中的高斯参量进行预测,实现GM-δ-GLMB的预测步骤;
(3)基于距离加权的量测合并策略。在预测与更新的迭代过程中假设轨迹的状态与量测关联关系未知的情况下,首先计算预测步骤2中的上一时刻的量测与当前时刻量测集中的统计距离,以此为依据通过查看χ2分布表设置阈值来提取与预测量测相关的量测集,然后,依据统计距离的大小对阈值内的量测赋予不同的权重,最后合并不同权重的量测,获取观测目标更新步骤对应的量测,从而基于距离加权的量测合并策略获取观测的地面目标更新步骤对应的多个传感器量测信息,用其更新地面多目标的伯努利参量。
(4)更新。在更新步骤中的δ-GLMB的空间分布密度函数仍然由高斯分量组成,并结合步骤3中合并后的传感器量测信息,利用CKF对δ-GLMB 的高斯分量进行更新;,得到本时刻最终输出的状态集合,对多个目标状态进行提取并估计目标运动轨迹,完成一次迭代;输出下一时刻多个传感器的量测合集,转至步骤S2进行下一次迭代,直至无关于下一刻的输入;
(5)高斯项的剪枝、合并。
(6)获取各个时刻的目标状态及目标运动轨迹,对多个目标状态进行提取并估计目标运动轨迹。
二、一个具体实施例如下:
本发明针对繁华闹市区,路况及地形复杂的情况,且存在大量的杂波及虚警,为了降低目标跟踪算法的计算复杂度且实时准确估计目标的运动状态及跟踪轨迹,本发明专利提出一种基于距离加权的量测组合策略的 CKF-δ-GLMB方法:获取目标的初始运动状态,基于步骤2至步骤4求得更新步骤的高斯参量。
步骤1:本发明采用GM来实现δ-GLMB,GM是将用于传递的伯努利参数利用多个高斯项加权求和的方式来代替,以解决预测与更新步骤中带有有限集积分运算而无法求得解析解的问题。我们第一步要定义并解释一些符号,按照以前的惯例,小写字母x可以表示单目标状态,例如:大写字母X可表示多目标状态,小写字母z和大写字母Z可以区分单、多目标的传感器量测值。先验概率密度中标签l对应的伯努利参数的GM形式为:
Figure RE-GDA0003988893070000131
其中,k为观测时间,
Figure RE-GDA0003988893070000132
表示先验概率密度中标签l对应参数所分解成的高斯项个数,
Figure RE-GDA0003988893070000133
表示k-1时刻的第j个目标的高斯项状态均值,j表示先验概率密度中标签l对应参数所分解成的高斯项个数,也就是第j个假设轨迹。
Figure RE-GDA0003988893070000134
表示其对应协方差,
Figure RE-GDA0003988893070000135
表示k-1时刻标签l第j个GM项的权重,在实际仿真中,此值为随机非负数,且满足
Figure RE-GDA0003988893070000136
表示k-1 时刻的第j个高斯项的状态均值为
Figure RE-GDA0003988893070000137
协方差为
Figure RE-GDA0003988893070000138
的高斯概率密度;y表示自变量参数。因此对高斯项的传递过程即为其参数从初始化到预测步骤最后到下一时刻的更新步骤的过程。
步骤2:为了提升GM-δ-GLMB在非线性系统中的滤波效果,本发明引入 CKF来对GM-δ-GLMB的高斯参量进行预测和更新,预测具体步骤如下:
(1)生成容积点集:
Figure RE-GDA0003988893070000139
其中,s=2n,n为系统状态维度,[1]=[In×n,In×n],I是n维单位矩阵,[1]t是单位矩阵的第t列向量。CKF的容积点集{ξt,wt},利用三阶球面径向规则计算标准高斯加权积分,即
Figure RE-GDA0003988893070000141
ξt与wt都是容积点集的参数,分别为基本容积点及权值,下标t代表单位矩阵的第t列向量。
(2)预测步骤,实现高斯参量
Figure RE-GDA0003988893070000142
Figure RE-GDA0003988893070000143
的传递:
Cholesky分解:
Figure RE-GDA0003988893070000144
Figure RE-GDA0003988893070000145
是对公式(1)组成的目标轨迹中标签l在k-1时刻的第j个高斯项的协方差矩阵进行Cholesky分解后的值。
产生容积点:
Figure RE-GDA0003988893070000146
Figure RE-GDA0003988893070000147
表示标签l在k-1时刻的第j个高斯项的容积点。
计算容积点经非线性状态方程传播后的值:
Figure RE-GDA0003988893070000148
估计状态预测值:
Figure RE-GDA0003988893070000149
计算预测误差协方差矩阵:
Figure RE-GDA00039888930700001410
其中,chol(·)表示对协方差矩阵进行cholesky分解,f表示状态转移矩阵,
Figure RE-GDA0003988893070000151
表示标签l在k-1时刻的第j个高斯项的过程噪声协方差,过程噪声是相互独立的,服从均值为0、协方差矩阵为R的高斯分布。
步骤3:在标准的δ-GLMB算法模型中,每个目标最多只能产生一个量测,而一个传感器量测也只能匹配一个真实目标,所以在GM-δ-GLMB算法中,属于同一个目标的假设轨迹最多只能由一个量测与其关联。因此,依据传感器量测集合建模的RFS形式,采用能够体现其统计特性的统计距离,基于最近邻法的思想,将与假设轨迹量测的预测量测值相近的量测作为可能与假设轨迹相关的量测。以此定义
Figure RE-GDA0003988893070000152
Figure RE-GDA0003988893070000153
的统计距离为:
Figure RE-GDA0003988893070000154
Figure RE-GDA0003988893070000155
为标签l在k-1时刻的第j个高斯项的新息协方差矩阵,公式(15) 可求得新息协方差矩阵。
根据上式(8),通过选定特定的关联阈值,第j个假设轨迹关联的量测可用式子(9)的统计门限提取:
Figure RE-GDA0003988893070000156
式中Th为第k时刻实现预定的阈值。
由上述式子求得第j个假设轨迹的相关量测,为了减小阈值内的杂波对状态更新的影响,依据阈值内量测统计距离的大小,合并每个量测分配相应的权值,以减小关联量测附近的杂波对目标状态估计的影响。
Figure RE-GDA0003988893070000161
其对应的权值可以描述为:
Figure RE-GDA0003988893070000162
Figure RE-GDA0003988893070000163
表示单目标在k时刻标签l第j个假设轨迹的相关量测,
Figure RE-GDA0003988893070000164
表示多目标在k时刻标签l第j个假设轨迹相关量测。
式中:
Figure RE-GDA0003988893070000165
为归一化因子,显然
Figure RE-GDA0003988893070000166
由上式可知,
Figure RE-GDA0003988893070000167
可以表示为:
Figure RE-GDA0003988893070000168
依据以上推导,可求得组合后的量测集
Figure RE-GDA0003988893070000169
并代入更新步骤中式(19),可求得假设轨迹的状态均值。
其中,
Figure RE-GDA00039888930700001610
表示k-1时刻单目标标签l第j个高斯项的传感器量测,
Figure RE-GDA00039888930700001611
表示k时刻单目标标签l第j个高斯项的传感器量测。
其中Th的设置可依据查看χ2分布表来获得,太小的Th会使得大量的量测不能用于更新操作,导致估计误差的增加;太大的Th则会使计算复杂度过大。
步骤4:CKF-GM-δ-GLMB的更新步骤,实现高斯参量从预测参数到下一时刻更新部分参数的传递:
更新步Cholesky分解:
Figure RE-GDA0003988893070000171
Figure RE-GDA0003988893070000172
表示对预测结果的第j个高斯项的协方差矩阵进行Cholesky分解值。
计算容积点:
Figure RE-GDA0003988893070000173
Figure RE-GDA0003988893070000174
为基本容积点。
计算容积点经传感器量测非线性函数传播后的值:
Figure RE-GDA0003988893070000175
其中,
Figure RE-GDA0003988893070000176
指标签l在k时刻的第j个高斯项的容积点经传感器量测非线性函数传播后的值。
计算新息协方差:
Figure RE-GDA0003988893070000177
其中,
Figure RE-GDA0003988893070000178
为标签l在k时刻的第j个高斯项的量测噪声协方差。
计算量测预测值:
Figure RE-GDA0003988893070000179
计算互协方差:
Figure RE-GDA00039888930700001710
计算卡尔曼增益:
Figure RE-GDA0003988893070000181
更新状态均值:
Figure RE-GDA0003988893070000182
更新协方差矩阵:
Figure RE-GDA0003988893070000183
其中,h(·)表示量测转移函数,
Figure RE-GDA0003988893070000184
是量测噪声协方差,量测噪声是相互独立的,服从均值为0、协方差矩阵为Q的高斯分布。由此,利用CKF实现了高斯参量预测步到更新步
Figure RE-GDA0003988893070000185
的转变。
输入下一时刻的多目标量测集合,转到步骤2进行迭代,直到无关于下一时刻的输入,整个算法结束。
步骤5:高斯项的剪枝、合并。舍去低于设定的生存概率的最低门限的伯努利分量,保留高于该门限的伯努利分量对应的高斯分量,对这些高斯分量进行剪枝,合并相距较近的高斯分量。
步骤6:对多个目标状态进行提取并估计目标运动轨迹。采用最大后验估计从势分布中估计多目标个数,从估计目标数相同的假设中选取权重最大假设的多目标均值与标签集合,分别作为估计目标状态及目标的估计轨迹。
三、本发明的效果可通过以下仿真进一步说明
本次仿真实验以地面多目标为研究对象,仿真实验中有八个地面机动目标在传感器监测区域[-2000m,2000m]×[0m,2000m]内运动。地面机动目标的运动状态为
Figure RE-GDA0003988893070000191
其中(xk,yk)为目标分别在x轴方向和y轴方向的位置信息,
Figure RE-GDA0003988893070000192
为分别沿着x轴方向和y轴方向的速度信息,(tstart,tstop)为目标的出现和消亡时间,wturn=2*p/180为转弯速率。
全程跟踪过程时长连续进行100s,采样时间间隔是T=1s,地面多目标存活概率为pS(x)=0.99,目标检测概率为pD,k(x)=0.98。目标1在1~100s做匀变速转弯运动,初始状态为X1=[1003.8676;-10;1488.2543;-10;wturn/8];目标2在10~100s做匀变速转弯运动,初始状态为X2=[-255.8857;20; 1011.4102;3;-wturn/3];目标3在10~66s做匀变速转弯运动,初始状态为X3=[-1500;43;250;0;0];目标4在20~80s做匀变速转弯运动,初始状态为X4=[246.1324;11;738.9253;5;wturn/4];目标5在40~100s做匀变速转弯运动,初始状态为X5=[-242.6194;-12;993.2007;-12;wturn/2];目标6在40~100s做匀变速转弯运动,初始状态为X6=[1000;0;1500;-10; wturn/4];目标7在40~80s做匀变速转弯运动,初始状态为X7=[250;-50; 750;0;-wturn/4];目标8在60~100s做匀变速转弯运动,初始状态为X8=[250;-40;750;25;wturn/4],八个地面目标的真实轨迹如图2所示,其中O与D分别代表每个地面目标的起始位置与终止位置。本发明算法在非线性模型系统下的目标运动模型:转弯模型为匀变速转弯运动CT,其状态转移方程为:
Figure RE-GDA0003988893070000201
其系统非线性量测方程为:
Figure RE-GDA0003988893070000202
其中,vk和wk分别为独立的过程噪声和观测噪声,它们遵循均值为0、协方差矩阵分别为Rk、Qk的高斯分布,过程噪声标准差σw=15m/s2,量测噪声协方差σθ=(π/180)rad/s,σr=5m,转速噪声标准差σu=π/180rad/s。
杂波模型服从强度为κk=λVu(z)的泊松RFS,其中V表示传感器监视区域的观测面积;u(z)表示该传感器监视区域的均匀概率分布函数,λ为杂波密度。仿真实验在MATLABR2016a环境下进行,处理器为Intel Core i7-7700,运行内存为8GB,本发明进行500次蒙特卡洛仿真实验来验证所提算法的有效性。在仿真场景中,新生目标模型服从多伯努利RFS,即
Figure RE-GDA0003988893070000203
其中,
Figure RE-GDA0003988893070000204
为目标新生概率,空间分布密度函数
Figure RE-GDA0003988893070000205
具体参数如下:
Figure RE-GDA0003988893070000206
Figure RE-GDA0003988893070000207
四、仿真实验结果分析
图3为本发明所提算法在杂波密度为λ=1.0×10-6m-2时的目标滤波结果图,
本发明所提算法在x轴和y轴方向上对地面机动目标的跟踪效果如图3所示。
从图3、图4中可看出本发明目标跟踪算法能够准确地估计目标的运动轨迹,对目标进行准确跟踪,并能够很好的区分不同目标轨迹。本发明算法仍选用最优子模式分配OSPA距离作为评价多目标跟踪算法性能的准则指标。本发明取p=1,c=100。
五、强杂波密度下改进算法的结果分析
在较偏僻的郊区场景及旷野环境下,跟踪地面目标的杂波密度比较小,提高跟踪算法的跟踪精度是重中之重。为了提高δ-广义标签多伯努利算法的跟踪精度,采用非线性滤波算法均方根容积卡尔曼滤波器(SCKF)与δ-GLMB 相结合。因多目标跟踪滤波算法δ-GLMB在跟踪滤波迭代过程中可能会随着运算次数的增加出现协方差矩阵非正定情况,这导致滤波器中不能使用均方根法,Cholesky分解,并导致滤波迭代结束,于是均方根容积卡尔曼滤波器 (SCKF)为了解决上述问题,以状态估计误差协方差的均方根传递预测与更新误差协方差,保持了协方差矩阵的正定性与对称性,提高δ-GLMB跟踪算法的稳定性,使用协方差矩阵的均方根作为递推多伯努利参数的之一。容积卡尔曼滤波器(CKF)虽然没有SCKF算法的稳定性及精确度高,但是CKF可以在满足基本精度保障的基础上,对非线性环境下的多目标进行跟踪。而且因为其采用三阶球面-径向容积准则近似估计非线性滤波算法中的概率积分,通过2n个容积点进行加权求和,在运算过程中没有使用协方差矩阵的均方根迭代运行,没有在预测更新环节的每个采样周期内进行QR分解来实现算法递归,因此其计算复杂度小于SCKF。因此本方法在面向强杂波密度场景下采用 CKF-δ-GLMB算法。
在强杂波环境下即在高杂波密度为λ=5.0×10-6m-2下进行与算法 SCKF-GM-δ-GLMB的OSPA距离对比实验仿真图如图5及势估计误差对比图如下图6。
由图5及图6算法对比图可以看出两种算法的目标状态信息估计情况相差不多,在势分布估计也就是目标数目估计情况两种算法的估计性能也相差无几,可以证明本发明算法与SCKF-GM-δ-GLMB算法的跟踪地面目标的性能效果相差不多。
表1不同杂波密度下运行时间比较
Figure RE-GDA0003988893070000221
实验结果分析:从表1可以看出,本发明算法CKF-GM-δ-GLMB的改进算法在不同高杂波密度下的运行时间都要明显低于SCKF-GM-δ-GLMB,对三个高杂波密度下的运算耗时综合比较,CKF-GM-δ-GLMB的改进算法平均耗时相较于SCKF-GM-δ-GLMB来说减少了20.92%。
本发明提出了一种改进CKF-GM-δ-GLMB多目标跟踪滤波算法。在强杂波密度场景下,本算法能够实现多目标实时跟踪,在计算效率方面有明显的提升,且从其运行时间随杂波密度增加的涨幅来看,本发明的运行时间涨幅略小,进一步验证了本发明的有效性。
本发明中,所谓的地面多目标跟踪(ground multi-target tracking),是指通过单个或者多个传感器获取地面多目标(如坦克、车辆等)的量测信息,通过一系列优化的综合滤波处理,实时采集、估计传感器监测区域内的地面目标数量、位置、速度和加速度等目标状态信息或运动轨迹,进而识别地面目标属性意图,分析和估计威胁评估的发展态势。
本发明中,所谓的随机有限集(random finite set,RFS),又称集值随机集,是一个以没有顺序的集合为元素的随机变量,也就是有限集值的随机变量。对于一个随机有限集,其数目与元素都随时间变化的,且取值本身随机无序,各元素间具有独立性,因此随机集中的元素具有很好的统计特性。
本发明中,所谓的容积卡尔曼滤波(Cubature Kalman filter,CKF) 是由加拿大学者Arasaratnam和Haykin在2009年首次在硕士学位论文提出,采用三阶球面-径向容积准则近似估计非线性滤波算法中的概率积分,通过2n个容积点进行加权求和,在运算过程中没有使用协方差矩阵的均方根迭代运行,没有在预测更新环节的每个采样周期内进行QR分解来实现算法递归,是理论上当前最接近贝叶斯滤波的近似算法,是解决非线性系统状态估计的强有力工具。
本发明中,所谓的δ-广义标签多目标多伯努利跟踪滤波 (δ-GeneralizedLabeled Multi-Bernoulli Filter,δ-GLMB),为 Ba-Tuon Vo等人将多假设跟踪思想(MHT)与随机有限集(RFS)理论相结合,提出了δ-广义标签多目标多伯努利滤波算法(δ-GLMB),该算法实现了真正意义上的目标跟踪,传递不同跟踪标签的多个假设集合。δ-GLMB算法无需共轭近似处理,具有更高的跟踪精度,抗干扰性能良好且能直接估计出目标的运动轨迹。
虽然本发明已利用上述较佳实施例进行说明,然其并非用以限定本发明的保护范围,任何本领域技术人员在不脱离本发明的精神和范围之内,相对上述实施例进行各种变动与修改仍属本发明所保护的范围,因此本发明的保护范围以权利要求书所界定的为准。

Claims (6)

1.地面多目标跟踪状态及轨迹确定方法,其特征在于包括如下步骤:
步骤S1:K=0时刻,初始化参数,获取地面多机动目标状态均值、多目标状态协方差矩阵、多目标非线性运动模型及杂波参数,得到初始分量集;采用GM来实现δ-GLMB,将用于传递的伯努利参数利用多个高斯项加权求和的方式来代替;
步骤S2:K大于或等于1时,通过非线性滤波器CKF对GM-δ-GLMB中预测空间分布密度函数的存活部分所分解的高斯项进行预测,完成对初始化的地面多目标状态集合信息的预测;
步骤S3:获取步骤S2中上一时刻的量测与当前时刻量测集中的统计距离,以此为依据通过查看χ2分布表设置阈值来提取与预测量测相关的量测集,然后,依据统计距离的大小对阈值内的量测赋予不同的权重,最后合并不同权重的量测,从而基于距离加权的量测合并策略获取观测的地面目标更新步骤对应的多个传感器量测信息;
步骤S4:结合步骤S3所获取的传感器量测信息,利用CKF对δ-GLMB的高斯分量进行更新,得到本时刻最终输出的状态集合,对多个目标状态进行提取并估计目标运动轨迹,完成一次迭代;输出下一时刻多个传感器的量测合集,转至步骤S2进行下一次迭代,直至无关于下一刻的输入;
步骤S5:对高斯项进行剪枝和合并;
步骤S6:获取各个时刻的目标状态及目标运动轨迹,对多个目标状态进行提取并估计目标运动轨迹。
2.如权利要求1所述的地面多目标跟踪状态及轨迹确定方法,其特征在于:
所述步骤S1中,采用GM来实现δ-GLMB时,先验概率密度中标签l对应的伯努利参数的GM形式为:
Figure FDA0003880542390000021
其中,k为观测时间,
Figure FDA0003880542390000022
表示先验概率密度中标签l对应参数所分解成的高斯项个数,j表示先验概率密度中标签l对应参数所分解成当前的高斯项,
Figure FDA0003880542390000023
表示k-1时刻的第j个目标的高斯项状态均值,
Figure FDA0003880542390000024
表示其对应协方差,
Figure FDA0003880542390000025
表示k-1时刻标签l第j个高斯项的权重,在实际仿真中,此值为随机非负数,且满足
Figure FDA0003880542390000026
表示k-1时刻的第j个高斯项的状态均值为
Figure FDA0003880542390000027
协方差为
Figure FDA0003880542390000028
的高斯概率密度;y表示自变量参数。
3.如权利要求2所述的地面多目标跟踪状态及轨迹确定方法,其特征在于,所述步骤S2包括如下步骤:
步骤S21:生成容积点集:
Figure FDA0003880542390000029
其中,s=2n,n为系统状态维度,[1]=[In×n-In×n],I是n维单位矩阵,[1]t是单位矩阵的第t列向量。CKF的容积点集{ξt,wt}},利用三阶球面径向规则计算标准高斯加权积分,即
Figure FDA0003880542390000031
ξt与wt都是容积点集的参数,分别为基本容积点及权值,下标t代表单位矩阵的第t列向量;
步骤S22:预测步骤,实现高斯参量
Figure FDA0003880542390000032
Figure FDA0003880542390000033
的传递,包括:
步骤S221,Cholesky分解:
Figure FDA0003880542390000034
Figure FDA0003880542390000035
是对公式(1)组成的目标轨迹中标签l在k-1时刻的第j个高斯项的协方差矩阵进行Cholesky分解后的值;
步骤S222,产生容积点:
Figure FDA0003880542390000036
Figure FDA0003880542390000037
表示标签l在k-1时刻的第j个高斯项的容积点;
步骤S223,计算容积点经非线性状态方程传播后的值:
Figure FDA0003880542390000038
步骤S224,估计状态预测值:
Figure FDA0003880542390000039
步骤S225,计算预测误差协方差矩阵:
Figure FDA00038805423900000310
其中,chol(·)表示对协方差矩阵进行cholesky分解,f表示状态转移矩阵,
Figure FDA0003880542390000041
表示标签l在k-1时刻的第j个高斯项的过程噪声协方差,过程噪声是相互独立的,服从均值为0、协方差矩阵为R的高斯分布。
4.如权利要求3所述的地面多目标跟踪状态及轨迹确定方法,其特征在于,所述步骤S3中,通过如下步骤推导出单目标的量测集:
步骤S31:定义
Figure FDA0003880542390000042
Figure FDA0003880542390000043
的统计距离为:
Figure FDA0003880542390000044
其中,
Figure FDA0003880542390000045
为标签l在k-1时刻的第j个高斯项的新息协方差矩阵;
步骤S32:通过选定特定的关联阈值,利用下式的统计门提取第j个假设轨迹关联的量测:
Figure FDA0003880542390000046
其中,Th为第k时刻实现预定的阈值;
步骤S33:
Figure FDA0003880542390000047
其对应的权值被描述为:
Figure FDA0003880542390000048
其中,
Figure FDA0003880542390000049
表示单目标在k时刻标签l第j个假设轨迹的相关量测,
Figure FDA00038805423900000410
表示多目标在k时刻标签l第j个假设轨迹相关量测;
Figure FDA00038805423900000411
为归一化因子;
步骤S34:将量测集表示为:
Figure FDA0003880542390000051
其中,
Figure FDA0003880542390000052
表示k-1时刻单目标标签l第j个高斯项的传感器量测,
Figure FDA0003880542390000053
表示k时刻单目标标签l第j个高斯项的传感器量测。
5.如权利要求4所述的地面多目标跟踪状态及轨迹确定方法,其特征在于,所述步骤S31中,式(8)中的新息协方差矩阵
Figure FDA0003880542390000054
通过如下式子得到:
步骤S311,更新步Cholesky分解:
Figure FDA0003880542390000055
其中,
Figure FDA0003880542390000056
表示对预测结果的标签l第j个高斯项的协方差矩阵进行Cholesky分解值;
步骤S312,计算容积点:
Figure FDA0003880542390000057
其中,
Figure FDA0003880542390000058
为标签l在k时刻的第j个高斯项的基本容积点;
步骤S313,计算容积点经传感器量测非线性函数传播后的值:
Figure FDA0003880542390000059
其中,
Figure FDA00038805423900000510
指标签l在k时刻的第j个高斯项的容积点经传感器量测非线性函数传播后的值;
步骤S314,计算新息协方差:
Figure FDA0003880542390000061
其中,
Figure FDA0003880542390000062
为标签l在k时刻的第j个高斯项的量测噪声协方差。
6.如权利要求5所述的地面多目标跟踪状态及轨迹确定方法,其特征在于,所述步骤S4包括:
步骤S41:更新步Cholesky分解,操作同于步骤S311;
步骤S42:计算容积点,操作同于步骤S312;
步骤S43:计算容积点经传感器量测非线性函数传播后的值,操作同于步骤S313;
步骤S44:计算信息协方差,操作同于步骤S314;
步骤S45,计算量测预测值:
Figure FDA0003880542390000063
步骤S46,计算互协方差:
Figure FDA0003880542390000064
步骤S47,计算卡尔曼增益:
Figure FDA0003880542390000065
步骤S48,更新状态均值:
Figure FDA0003880542390000066
步骤S49,更新协方差矩阵:
Figure FDA0003880542390000071
至此,实现高斯参量预测步到更新步
Figure FDA0003880542390000072
的转变;
其中,h(·)表示量测转移函数,
Figure FDA0003880542390000073
是量测噪声协方差,量测噪声是相互独立的,服从均值为0、协方差矩阵为Q的高斯分布;
输出下一时刻的多目标量测集合,转到步骤S2进行迭代,直到无关于下一时刻的输入时,转至步骤S5。
CN202211227444.XA 2022-10-09 2022-10-09 地面多目标跟踪状态及轨迹确定方法 Pending CN115619825A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211227444.XA CN115619825A (zh) 2022-10-09 2022-10-09 地面多目标跟踪状态及轨迹确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211227444.XA CN115619825A (zh) 2022-10-09 2022-10-09 地面多目标跟踪状态及轨迹确定方法

Publications (1)

Publication Number Publication Date
CN115619825A true CN115619825A (zh) 2023-01-17

Family

ID=84861423

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211227444.XA Pending CN115619825A (zh) 2022-10-09 2022-10-09 地面多目标跟踪状态及轨迹确定方法

Country Status (1)

Country Link
CN (1) CN115619825A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117724087A (zh) * 2024-02-07 2024-03-19 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117724087A (zh) * 2024-02-07 2024-03-19 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法
CN117724087B (zh) * 2024-02-07 2024-05-28 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Similar Documents

Publication Publication Date Title
Li et al. Joint data association, registration, and fusion using EM-KF
Svensson et al. Set JPDA filter for multitarget tracking
CN107831490A (zh) 一种改进的多扩展目标跟踪方法
CN111127523B (zh) 基于量测迭代更新的多传感器gmphd自适应融合方法
Beard et al. A generalised labelled multi-Bernoulli filter for extended multi-target tracking
Dong et al. Maneuvering multi-target tracking based on variable structure multiple model GMCPHD filter
Clark et al. Data association for the PHD filter
Yang et al. Linear-time joint probabilistic data association for multiple extended object tracking
Lan et al. Joint target detection and tracking in multipath environment: A variational Bayesian approach
Liangqun et al. Maximum entropy fuzzy clustering with application to real-time target tracking
CN113673565B (zh) 多传感器gm-phd自适应序贯融合多目标跟踪方法
CN111562571A (zh) 一种未知新生强度的机动多目标跟踪与航迹维持方法
CN109214432B (zh) 一种多传感器多目标联合检测、跟踪与分类方法
Yang et al. A novel track maintenance algorithm for PHD/CPHD filter
Qin et al. Application of an efficient graph-based partitioning algorithm for extended target tracking using GM-PHD filter
CN115619825A (zh) 地面多目标跟踪状态及轨迹确定方法
CN116630370A (zh) 多模型pbp-tpmb机动扩展目标跟踪方法
Li et al. Labeled multi-Bernoulli filter based multiple resolvable group targets tracking with leader–follower model
Sunyoung et al. Cardinality compensation method based on information-weighted consensus filter using data clustering for multi-target tracking
CN111340853B (zh) 基于ospa迭代的多传感器gmphd自适应融合方法
Nandakumaran et al. Improved multitarget tracking using probability hypothesis density smoothing
CN113537302B (zh) 一种多传感器随机化数据关联融合方法
Vu et al. Particle Markov chain Monte Carlo for Bayesian multi-target tracking
García-Fernández et al. Continuous-discrete trajectory PHD and CPHD filters
CN113486300A (zh) 一种无人驾驶车辆多目标跟踪方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination