CN111562571B - 一种未知新生强度的机动多目标跟踪与航迹维持方法 - Google Patents

一种未知新生强度的机动多目标跟踪与航迹维持方法 Download PDF

Info

Publication number
CN111562571B
CN111562571B CN202010466368.2A CN202010466368A CN111562571B CN 111562571 B CN111562571 B CN 111562571B CN 202010466368 A CN202010466368 A CN 202010466368A CN 111562571 B CN111562571 B CN 111562571B
Authority
CN
China
Prior art keywords
target
particles
measurement
particle
time
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.)
Active
Application number
CN202010466368.2A
Other languages
English (en)
Other versions
CN111562571A (zh
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.)
Jiangnan University
Original Assignee
Jiangnan University
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 Jiangnan University filed Critical Jiangnan University
Priority to CN202010466368.2A priority Critical patent/CN111562571B/zh
Publication of CN111562571A publication Critical patent/CN111562571A/zh
Application granted granted Critical
Publication of CN111562571B publication Critical patent/CN111562571B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/66—Radar-tracking systems; Analogous systems
    • G01S13/70—Radar-tracking systems; Analogous systems for range tracking only

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种未知新生强度的机动多目标跟踪与航迹维持方法，属于智能信息处理技术领域。本发明方法在CPHD滤波框架下，引入参数自适应估计和粒子标识航迹关联技术以及新生目标识别策略，提出一种基于参数自适应CPHD滤波方法，以解决对复杂环境下新生目标强度未知，数目未知且时变的机动多目标跟踪的问题。本发明方法中将目标状态和时变的模型参数进行联合在线估计，采用包含不同模型参数的粒子对系统模型进行融合估计，以提高对机动目标的适应能力；滤波过程中对所有粒子进行身份标识，实现了对于新生目标可以通过量测自动识别和对多目标的航迹管理；具有较强的鲁棒性和抗干扰能力，可以满足实际工程系统的设计需求，具有良好的工程应用价值。

Description

一种未知新生强度的机动多目标跟踪与航迹维持方法
技术领域
本发明属于智能信息处理技术领域，具体涉及一种未知新生强度的机动多目标跟踪与航迹维持方法。
背景技术
在多目标跟踪领域中，早期主要采用数据关联技术实现对多目标的跟踪，如联合概率数据关联(JPDA)和多假设跟踪(MHT)等，虽然这些方法在多目标跟踪中具有一定的效果，但由于存在复杂的数据关联运算，尤其是计算复杂度随着目标个数的增加而呈现指数增长，影响算法的实时性。此外，对数目未知且变化的多目标跟踪，存在目标数目及状态估计不准确的问题。
近年来，随机有限集(Random Finite Set，RFS)理论在对数目未知且变化的多目标跟踪中取得了一定优势，分别对目标状态和观测进行随机集建模，可避免复杂的数据关联运算。自Mahler教授提出概率假设密度(Probability hypothesis density，PHD)滤波之后，随机有限集理论在目标跟踪领域得到了广泛应用，随后，他放宽目标数目服从泊松分布的限制，提出势均衡的PHD(CPHD)滤波方法，该方法能够联合估计目标随机集的概率假设密度和目标数目的概率分布，相比PHD滤波方法，进一步提高了对数目未知且变化的多目标跟踪性能。典型的闭合求解方式有：高斯混合PHD/CPHD和粒子PHD/CPHD，但这些方法难以实现对复杂环境下，新生强度未知，量测噪声未知且任意机动的多目标跟踪。
针对机动目标跟踪，主要包含基于单模型和多模型的滤波方法，其中，单模型方法如Singer模型，当前统计模型和输入估计模型等，对机动目标的跟踪性能往往取决于对目标机动参数的选取，如机动频率或最大加速度等，如果参数选取不合适，将严重影响算法的跟踪性能；而多模型方法，则采用多个模型进行匹配滤波，对机动目标的跟踪性能通常取决于模型集的设定，如果模型集中的模型设置不合理，也将直接影响算法的跟踪性能，且模型集的数目也将直接影响算法的运算效率。
发明内容
为了克服现有方法存在的不足，本发明拟在CPHD滤波框架下，引入自适应参数估计(Adaptive parameter estimate,APE)以及新生目标识别策略，对未知的噪声参数和机动参数进行自适应估计，并采用粒子标记技术对每个目标进行身份标识，实现对复杂环境下任意机动多目标自适应跟踪和航迹管理。自适应参数估计是Christopher Nemeth等于2014年在《Sequential Monte Carlo Methods for State and Parameter Estimation inAbruptly Changing Environments》中提出，借助Liu and West(LW)滤波实现对任意时变参数估计，采用逆伽马(Inverse-Gamma，IG)分布近似估计未知的静态参数，本发明中引入该方法来估计机动多目标跟踪场景中的时变机动参数和未知的量测噪声。此外，滤波过程中对所有粒子进行身份标识，并通过相似度来实现了对新生目标的识别，实现了对多目标的航迹管理。
本发明在CPHD滤波框架下，引入参数自适应估计、新生目标识别技术和粒子标识航迹关联技术，提出一种基于参数自适应CPHD滤波方法，以解决对复杂环境下数目未知且时变、新生强度未知的机动多目标跟踪。发明方法中将目标状态和时变的模型参数进行联合在线估计，采用包含不同模型参数的粒子对系统模型进行融合估计，以提高对机动目标的适应能力；此外，滤波过程中对所有粒子进行身份标识，并通过相似度来实现了对新生目标的识别，实现了对多目标的航迹管理。
本发明涉及一种用于数目未知且变化的机动多目标跟踪及航迹管理方法。具体地说是一种基于参数自适应粒子势概率假设密度(CPHD)滤波方法，在CPHD的滤波框架下，对未知的机动参数进行估计，实现对复杂环境下机动多目标跟踪。该方法可广泛用于雷达目标探测与跟踪、武器精确制导、低空突防，以及人机交互、智能交通管制、无人驾驶等应用领域。
实现本发明的关键技术是：将自适应参数估计技术融入到CPHD滤波框架中，以解决任意机动多目标跟踪的问题，其中，系统模型中包括未知的静态参数和时变参数，如量测噪声方差和目标时变的机动参数等，本发明中采用逆伽马分布近似量测噪声的后验分布，采用自适应LW滤波估计时变的机动参数，并采用新生目标识别技术和粒子标记技术实现对未知新生目标强度的多目标的航迹跟踪。
本发明的第一个目的是提供一种未知新生强度的多目标跟踪方法，所述方法是将采集的目标粒子按照目标物进行分类，获得粒子集
Figure GDA0003503971590000021
其中，k表示当前时刻，
Figure GDA0003503971590000022
表示分类数目；选择
Figure GDA0003503971590000023
中同类的所有粒子，计算每一类粒子的位置均值
Figure GDA0003503971590000024
计算每一类粒子位置均值与当前时刻每个量测位置的距离
Figure GDA0003503971590000025
其中
Figure GDA0003503971590000026
Mk表示当前时刻量测的个数。根据距离最短的
Figure GDA0003503971590000027
进行判断存活目标判断，表明第i类目标粒子与当前的第j个量测对应，k时刻未配对的量测记为Λk；计算每一类粒子位置均值与k+1时刻每个量测位置的距离
Figure GDA0003503971590000028
Mk+1表示k+1时刻量测的个数，根据距离最短的
Figure GDA0003503971590000029
进行判断存活目标判断，如果距离
Figure GDA0003503971590000031
小于阈值U，表明第i类目标粒子与当前的第j个量测对应，k+1时刻未配对的量测记为Λk+1，对于k时刻剩下的量测，计算与k+1时刻剩下量测位置的距离
Figure GDA0003503971590000032
其中，
Figure GDA0003503971590000033
表示k时刻量测集中没有和目标匹配的量测个数，
Figure GDA0003503971590000034
表示k+1时刻量测集中没有和目标匹配的量测的个数，若距离
Figure GDA0003503971590000035
小于阈值U，则k时刻第i个未关联的量测为新生目标，并根据该量测采样N个新生目标粒子，k时刻如果没有被关联的量测，判断为杂波。
在本发明的一种实施方式中，所述的阈值U是指根据多个目标的最大运动速度和采样间隔的乘积设定最大距离，反映当前时刻量测与下一时刻量测关联情况下的最大距离，当小于最大距离代表有新生目标生成。
本发明的第二个目的是提供一种机动多目标跟踪或航迹管理方法，所述方法包括如下步骤：
(1)初始化步骤：
(1a)初始时刻，假设有n0个目标，每个目标采样N个粒子，总粒子数为L0＝N×n0。初始目标状态集为X0，目标存在概率为Ps，目标检测概率为Pd，势分布为p0；
(1b)初始化采样粒子并对粒子进行标识，设P0(X0,θ0,ξ0)为先验的联合概率分布。从P0(X0,θ0,ξ0)中采用初始粒子集
Figure GDA0003503971590000036
每个粒子的权重设为
Figure GDA0003503971590000037
初始时刻，粒子标记表示为
Figure GDA0003503971590000038
初始时刻粒子的分类可表示为
Figure GDA0003503971590000039
Figure GDA00035039715900000310
表示分类数目；
(2)当k≥1时，产生2Lk-1个粒子，其中，从后验分布
Figure GDA00035039715900000311
中产生Lk-1个预测粒子
Figure GDA00035039715900000312
权重为
Figure GDA00035039715900000313
其中，i＝1,2,…,Lk-1，
Figure GDA00035039715900000314
是时变参数粒子
Figure GDA00035039715900000315
的均值，Zk是量测集；然后从从后验分布
Figure GDA00035039715900000316
中产生其余的Lk-1个粒子，权重为
Figure GDA00035039715900000317
其中，
Figure GDA00035039715900000318
为是从时变参数服从的先验分布Pξ(ξ0)中采样获得，i＝Lk-1+1,Lk-1+2,…,2Lk-1；
(3)计算量测噪声协方差：
(3a)假设
Figure GDA00035039715900000319
为静态参数粒子
Figure GDA00035039715900000320
的充分统计量，参数粒子通过
Figure GDA00035039715900000321
获得，其中～表示采样，P()表示分布，
Figure GDA00035039715900000322
表示粒子
Figure GDA00035039715900000323
从分布
Figure GDA00035039715900000324
中采样获得。考虑静态未知参数Vk为未知的量测噪声方差，其共轭先验可采用参数为a和b的逆伽马分布IG(a,b)近似。设定逆伽马分布参数，
Figure GDA0003503971590000041
Figure GDA0003503971590000042
其中，l＝1,…,d，d为量测噪声协方差R的维度；
(3b)计算量测噪声协方差
Figure GDA0003503971590000043
其中，n＝1,…,Lk-1为迭代次序，Lk-1为最大迭代次数；
(3c)如果n≤N，更新参数
Figure GDA0003503971590000044
返回步骤(3b)；
(4)计算每个粒子量测噪声协方差
Figure GDA0003503971590000045
利用量测噪声协方差计算前Lk-1个粒子对应于时变参数无变化点
Figure GDA0003503971590000046
预测似然成比例的另一个权重
Figure GDA0003503971590000047
其中i＝1,2,…,Lk-1；计算剩余Lk-1个粒子对应于时变参数变化点
Figure GDA0003503971590000048
的预测似然成比例的另一个权重
Figure GDA0003503971590000049
其中i＝Lk-1+1,Lk-1+2,…,2Lk-1；
(5)选择
Figure GDA00035039715900000410
中同类的所有粒子，计算每一类粒子的位置均值
Figure GDA00035039715900000411
计算每一类粒子位置均值与当前时刻每个量测位置的距离
Figure GDA00035039715900000412
其中
Figure GDA00035039715900000413
Mk表示当前时刻量测的个数；根据距离最短的
Figure GDA00035039715900000414
进行判断存活目标判断，表明第i类目标粒子与当前的第j个量测对应，k时刻未配对的量测记为Λk。计算每一类粒子位置均值与k+1时刻每个量测位置的距离
Figure GDA00035039715900000415
Mk+1表示k+1时刻量测的个数，根据距离最短的
Figure GDA00035039715900000416
进行判断存活目标判断，如果距离
Figure GDA00035039715900000417
小于阈值U，表明第i类目标粒子与当前的第j个量测对应，k+1时刻未配对的量测记为Λk+1，对于k时刻剩下的量测，计算与k+1时刻剩下量测位置的距离
Figure GDA00035039715900000418
其中，
Figure GDA00035039715900000419
表示k时刻量测集中没有和目标匹配的量测个数，
Figure GDA00035039715900000420
表示k+1时刻量测集中没有和目标匹配的量测的个数，若距离
Figure GDA00035039715900000421
小于阈值U，则k时刻第i个未关联的量测为新生目标，则根据当前对应量测位置采样N个新生目标粒子；
(6)选择2Lk-1个粒子中的Lk-1个粒子，将它们的索引表示为：l(i)∈{1,2,…,2Lk-1}，其中，i＝1,2,…,Lk-1。对于i＝1,2,…,Lk-1,选择索引li来自{1,2,…,Lk-1}的概率为
Figure GDA00035039715900000422
来自{Lk-1+1,Lk-1+2,…,2Lk-1}的概率为
Figure GDA00035039715900000423
其中β是发生突变的概率，并假设已知；
(7)预测势分布pk|k-1；
(8)预测粒子标记
Figure GDA00035039715900000424
和粒子分类
Figure GDA00035039715900000425
(9)更新目标权重
Figure GDA00035039715900000426
得到
Figure GDA00035039715900000427
(10)更新势分布pk；
(11)估计目标强度函数Dk；
(12)更新粒子标记
Figure GDA0003503971590000051
和粒子分类
Figure GDA0003503971590000052
(13)计算当前时刻总的目标数目nk；
(14)聚类获得同一目标簇的粒子，计算时变参数粒子
Figure GDA0003503971590000053
的均值
Figure GDA0003503971590000054
和协方差Vk-1的估计，更新时变参数
Figure GDA0003503971590000055
的值，得到
Figure GDA0003503971590000056
(15)重采样：
(15a)通过权重对粒子集
Figure GDA0003503971590000057
进行重采样得到新的粒子集
Figure GDA0003503971590000058
重采样之后每个粒子被赋予相同的权重
Figure GDA0003503971590000059
其中，Lk＝Lk-1+Jk,Jk是k时刻新生目标粒子数；
(15b)重采样后的粒子与它们的父辈粒子具有相同的标记,对重采样后的粒子进行标记，得到
Figure GDA00035039715900000510
其中，j＝1,2,…,Lk-1+Jk；
(16)通过聚类粒子获得目标状态，并且聚类中心是估计的状态
Figure GDA00035039715900000511
给新的分类分配新的标记
Figure GDA00035039715900000512
得到估计后的分类为
Figure GDA00035039715900000513
(17)通过粒子标记进行航迹关联得到每个目标的航迹，若下一时刻观测信息到达，转到步骤(2)进行迭代；否则，目标跟踪过程结束。
本发明具有以下优点：
(1)本发明中采用了带势分布的PHD(CPHD)滤波技术，有效实现了对数目未知且变化的视频多目标跟踪；且由于采用新生目标识别和粒子标记技术，可有效地识别出新生目标，且不需要新生强度等先验信息；
(2)本发明由于引入了自适应参数估计，对目标作机动运动，运动发生突变时能够正确的跟上目标，且可以有效解决由于机动运动而造成失跟的问题。
附图说明
图1是本发明的整体流程图；
图2是本发明方法中的机动多目标二维运动轨迹；
图3是本发明方法与PF-APE-PHD,PF-MM-PHD,PF-MM-CPHD算法平均OSPA距离对比图；
图4是本发明方法与PF-APE-PHD,PF-MM-PHD,PF-MM-CPHD算法平均目标数估计对比图；
图5是本发明方法具有不同的测量噪声方差平均OSPA距离对比图；
图6是本发明方法具有不同的测量噪声方差平均目标数估计对比图；
图7是本发明方法跟踪机动多目标的航迹图。
具体实施方式
为了更好的理解下述技术方案，作如下基础理论介绍：
1.CPHD滤波原理
对于多目标跟踪，当目标数未知或随时间变化时，目标数是一个离散随机变量，状态空间的维数也会随目标数的变化而变化。类似地，量测数也是一个随时间变化的离散随机变量，多目标的状态及量测数据集可以建模为随机有限集的形式，即：
Figure GDA0003503971590000061
Figure GDA0003503971590000062
其中，Xk为k时刻的目标状态集，Zk为量测集，F(X)、F(Z)分别是X和Z上的所有有限子集的集合，Nk为目标数，Mk为量测数，其中某些量测可能源于杂波。
若k-1时刻目标状态随机集为Xk-1，则k时刻的目标状态随机集Xk可表示为：
Figure GDA0003503971590000063
其中，Sk|k-1表示从k-1时刻到k时刻仍然存在的目标状态随机集，Bk|k-1表示k时刻衍生出来的目标状态随机集，Γk表示新生目标状态随机集。
目标量测随机集Zk可表示为
Figure GDA0003503971590000064
其中，Kk表示k时刻源于杂波的量测随机集，Θk表示源于真实目标的量测随机集。
令Dk|k-1和pk|k-1分别表示k-1时刻预测的多目标强度函数和势分布，Dk和pk表示k时刻多目标后验强度函数和势分布，CPHD滤波主要包含预测和更新步骤。
预测：
Figure GDA0003503971590000065
Dk|k-1(x)＝∫Psfk|k-1(x|x′)Dk-1(x′)dx′+∫βk|k-1(x|x′)Dk-1(x′)d(x′)+γk(x)
Figure GDA0003503971590000067
其中，pΓ,k是新生目标随机集的势分布，<.,.>表示内积运算，
Figure GDA0003503971590000071
表示二项式系数。
更新：
Figure GDA0003503971590000072
Figure GDA0003503971590000073
其中，
Figure GDA0003503971590000074
其中，
Figure GDA0003503971590000075
表示排列系数，|Z|表示量测集的势，pK,k是杂波随机集的势分布，Zk\{z}表示在Zk中除了z的剩余量测。
Ξk(D,Z)＝{<D,Ψk,z(x)>:z∈Z}
Figure GDA0003503971590000076
其中，ej是一个j阶的初等对称函数，即：
ej({ρ1,ρ2,…,ρm})＝(-1)jαm-j/αm
{ρ1,ρ2,…,ρm}是多项式αmxm+αm-1xm-1+...+α1x+α0的不同根。
2.自适应参数估计理论
Liu和West(LW)滤波器可用于联合识别静态参数和目标状态，它使用多元高斯分布的混合来近似并传播未知参数的边缘后验分布，将粒子引入到该滤波器中，获得的APE滤波器可以估计静态和时变参数。
APE方法中将参数联合后验分布P(xk,Φk|Z1:k-1)进行分解，即：
P(xk,Φk|Z1:k-1)＝P(xk|Z1:k-1,Φk)P(Φk|Z1:k-1)
其中，Φk＝[θk,ξk]，θk和ξk分别表示静态参数和时变参数向量，Φk的边缘预测分布可进一步表示为:
P(Φk|Z1:k-1)＝P(θk,ξk|Z1:k-1)＝P(θk|Z1:k-1,ξk)P(ξk|Z1:k-1)
P(θk|Z1:k-1,ξk)表示静态参数向量θk的预测分布，P(ξk|Z1:k-1)表示时变参数向量ξk的预测分布，可通过下式来近似计算：
Figure GDA0003503971590000081
其中，
Figure GDA0003503971590000082
表示均值为
Figure GDA0003503971590000083
协方差为h2Vk-1的高斯分布，h∈(0,1)表示缩放参数，
Figure GDA0003503971590000084
是第i个分量的权重，β为ξk在时间k处突然发生变化的概率。假设时变参量ξk在两个相邻变化点之间是分段常数，根据上式的定义，如果ξk没有突然变化，则其预测分布遵循N个分量的高斯混合模型，通过下式获得每个分量的均值和方差：
Figure GDA0003503971590000085
Figure GDA0003503971590000086
其中，
Figure GDA0003503971590000087
是ξk-1在k-1时刻估计的最小均方根误差，
Figure GDA0003503971590000088
表示时变参数向量ξk-1的第i个高斯分量。
Figure GDA0003503971590000089
是缩放因子，标准核平滑要求核分量以均值向量
Figure GDA00035039715900000810
为中心，这会导致后验分布过度分散，即高斯混合的协方差大于Vk-1，引入缩放因子可以迫使粒子
Figure GDA00035039715900000811
更接近于样本均值
Figure GDA00035039715900000812
从而保持相同的协方差Vk-1。k时刻参数突变时，时变参量ξk的预测分布将由初始分布Pξ(ξ0)决定。APE滤波中，可采用粒子滤波近似参数和状态联合的后验分布P(xk,Φk|Z1:k)，即：
Figure GDA00035039715900000813
假设在时间k-1时刻，后验分布由权重为
Figure GDA00035039715900000814
的N个粒子
Figure GDA00035039715900000815
近似表示，则在k时刻，每个粒子将根据
Figure GDA00035039715900000816
服从的两种不同分布而产生两个权值，即：
Figure GDA00035039715900000817
其中，
Figure GDA00035039715900000818
Figure GDA00035039715900000819
Figure GDA00035039715900000820
其中，
Figure GDA00035039715900000821
Figure GDA00035039715900000822
也即产生2N个粒子，
Figure GDA0003503971590000091
和
Figure GDA0003503971590000092
分别对应参数没有发生突变和发生突变时刻粒子的权值。最终从2N个粒子中分别按概率
Figure GDA0003503971590000093
和
Figure GDA0003503971590000094
进行粒子抽取，抽取N个粒子来近似表示参数和状态联合的后验分布P(xk,Φk|Z1:k)。
实施例1利用未知新生强度-参数自适应粒子势概率假设密度滤波方法跟踪机动多目标：
1.自适应参数CPHD滤波迭
为了简化表示，假设存活概率和检测概率都独立于目标状态向量和未知参数向量Φk，分别用PS和Pd表示，k-1时刻的联合后验概率假设密度表示为Dk-1(x,Φ)，则预测概率假设密度Dk|k-1(x,Φ)可表示为：
Figure GDA0003503971590000095
当k时刻量测已知时，联合后验概率假设密度可更新为：
Figure GDA0003503971590000096
其中，Φ是未知的参数向量，则Dk|k-1(x,Φ)很难直接计算获得，本方法采用粒子滤波技术获得其近似解。
2.具体实施步骤
参照图1，本发明的具体实施步骤包括如下：
步骤1.初始化步骤：
(1.1)初始时刻，假设有n0个目标，每个目标采样N个粒子，总粒子数为L0＝N×n0。初始目标状态集为X0，设P0(X0,θ0,ξ0)为先验的联合概率分布，从P0(X0,θ0,ξ0)中采用初始粒子集
Figure GDA0003503971590000097
每个粒子的权重设为
Figure GDA0003503971590000098
目标存在概率为Ps，目标检测概率为Pd，势分布为p0。粒子标记表示为
Figure GDA0003503971590000101
粒子的分类可表示为
Figure GDA0003503971590000102
Figure GDA0003503971590000103
表示分类数目。
步骤2.当k≥1时，目标预测：
(2.1)APE滤波过程中产生2Lk-1个粒子。其中，静态参数情况下从建议分布为
Figure GDA0003503971590000104
Figure GDA0003503971590000105
中采样Lk-1个粒子，即：
Figure GDA0003503971590000106
权重为：
Figure GDA0003503971590000107
其中，i＝1,2,…,Lk-1，
Figure GDA0003503971590000108
是转移概率密度函数。
另外的Lk-1个粒子是在未知时变参数情况下采样的粒子，假设时变参数服从先验先验分布Pξ(ξ0)，则未知时变参数可从先验分布中采用，即：
Figure GDA0003503971590000109
其中，i＝Lk-1+1,Lk-1+2,…,2Lk-1。则另外的Lk-1粒子可从建议分布
Figure GDA00035039715900001010
中采样获得，即：
Figure GDA00035039715900001011
权重为：
Figure GDA00035039715900001012
(2.2)假设
Figure GDA00035039715900001013
为静态参数粒子
Figure GDA00035039715900001014
的充分统计量，参数粒子通过
Figure GDA00035039715900001015
获得，其中～表示采样，P()表示分布，
Figure GDA00035039715900001016
表示粒子
Figure GDA00035039715900001017
从分布
Figure GDA00035039715900001018
中采样获得。考虑静态未知参数Vk为未知的量测噪声方差，其共轭先验可采用参数为a和b的逆伽马分布IG(a,b)近似。更新逆伽马分布参数，
Figure GDA00035039715900001019
Figure GDA00035039715900001020
其中，i＝1,2,…,Lk-1，l＝1,…,d，d为量测噪声协方差R的维度；
(2.3)计算量测噪声协方差
Figure GDA00035039715900001021
Figure GDA00035039715900001022
为迭代次序，N为最大迭代次数；
(2.4)如果n≤N，更新参数
Figure GDA00035039715900001023
Figure GDA00035039715900001024
返回步骤(2.3)；
(2.5)计算参数
Figure GDA0003503971590000111
Figure GDA0003503971590000112
计算每个粒子量测噪声协方差
Figure GDA0003503971590000113
利用量测噪声方差计算前Lk-1个粒子对应于时变参数无变化点
Figure GDA0003503971590000114
预测似然成比例的另一个权重
Figure GDA0003503971590000115
Figure GDA0003503971590000116
其中，i＝1,2,…,Lk-1，计算剩余Lk-1个粒子对应于时变参数变化点
Figure GDA0003503971590000117
的预测似然成比例的另一个权重
Figure GDA0003503971590000118
Figure GDA0003503971590000119
其中，i＝Lk-1+1,Lk-1+2,…,2Lk-1。
(2.6)从2Lk-1个粒子中选取Lk-1个，将它们的索引表示为：l(i)∈{1,2,…,2Lk-1}，i＝1,2,…,Lk-1。对于i＝1,2,…,Lk-1,选择索引li来自{1,2,…,Lk-1}的概率为
Figure GDA00035039715900001110
来自{Lk-1+1,…,2Lk-1}的概率为
Figure GDA00035039715900001111
其中β是发生突变的概率，并假设已知。如果li∈{1,2,…,Lk-1}，然后使用下式选择时变参数粒子：
Figure GDA00035039715900001112
使用下式：
Figure GDA00035039715900001113
Figure GDA00035039715900001114
将参数粒子复合成：
Figure GDA00035039715900001115
如果li∈{Lk-1+1,Lk-12,…,2Lk-1}，然后把时变参数粒子设置为：
Figure GDA00035039715900001116
使用下式：
Figure GDA00035039715900001117
Figure GDA00035039715900001118
将参数粒子复合成：
Figure GDA00035039715900001119
使用索引i＝1,2,…,Lk-1重新标识选择的粒子：
Figure GDA00035039715900001120
Figure GDA0003503971590000121
选择
Figure GDA0003503971590000122
中同类的所有粒子，计算每一类粒子的位置均值
Figure GDA0003503971590000123
其中，
Figure GDA0003503971590000124
其中，
Figure GDA0003503971590000125
计算每一类粒子位置均值与当前时刻每个量测位置的距离
Figure GDA0003503971590000126
其中
Figure GDA0003503971590000127
Mk表示当前时刻量测的个数。根据距离最短的
Figure GDA0003503971590000128
进行判断存活目标判断，表明第i类目标粒子与当前的第j个量测对应，k时刻未配对的量测记为Λk。计算每一类粒子位置均值与k+1时刻每个量测位置的距离
Figure GDA0003503971590000129
Mk+1表示k+1时刻量测的个数，根据距离最短的
Figure GDA00035039715900001210
进行判断存活目标判断，如果距离
Figure GDA00035039715900001211
小于阈值U，表明第i类目标粒子与当前的第j个量测对应，k+1时刻未配对的量测记为Λk+1，对于k时刻剩下的量测，计算与k+1时刻剩下量测位置的距离
Figure GDA00035039715900001212
其中，
Figure GDA00035039715900001213
表示k时刻量测集中没有和目标匹配的量测个数，
Figure GDA00035039715900001214
表示k+1时刻量测集中没有和目标匹配的量测的个数，若距离
Figure GDA00035039715900001215
小于阈值U，则k时刻第i个未关联的量测为新生目标，则根据当前对应量测位置采样N个新生目标粒子
Figure GDA00035039715900001216
总的新生粒子数为Jk。其中，i＝Lk-1+1,Lk-1+2,…Lk-1+Jk，阈值U为当前时刻量测与下一时刻量测关联情况下的最大距离。
(2.7)从2Lk-1个粒子中选取Lk-1个，将它们的索引表示为：l(i)∈{1,2,…2Lk-1}，i＝1,2,…,Lk-1。对于i＝1,2,…,Lk-1,选择索引li来自{1,2,…,Lk-1}的概率为
Figure GDA00035039715900001217
来自{Lk-1+1,…,2Lk-1}的概率为
Figure GDA00035039715900001218
其中β是发生突变的概率，并假设已知。如果li∈{1,2,…,Lk-1}，然后使用下式选择时变参数粒子：
Figure GDA00035039715900001219
使用下式：
Figure GDA00035039715900001220
Figure GDA00035039715900001221
将参数粒子复合成：
Figure GDA00035039715900001222
如果li∈{Lk-1+1,Lk-12,…,2Lk-1}，然后把时变参数粒子设置为：
Figure GDA0003503971590000131
使用下式：
Figure GDA0003503971590000132
Figure GDA0003503971590000133
将参数粒子复合成：
Figure GDA0003503971590000134
使用索引i＝1,2,…Lk-1重新标识选择的粒子：
Figure GDA0003503971590000135
Figure GDA0003503971590000136
(2.8)计算预测势分布pk|k-1：
Figure GDA0003503971590000137
Figure GDA0003503971590000138
其中，pΓ,k是新生目标随机集的势分布，<.,.>表示内积运算。
Figure GDA0003503971590000139
表示二项式系数。
(2.9)预测粒子标识：
Lk-1个预测粒子标记可表示为：
Figure GDA00035039715900001310
其中i＝1,2,…,Lk-1,预测阶段粒子的分类可表示为：
Figure GDA00035039715900001311
对于在预测阶段新生的Jk个粒子，给予新的粒子标记，即：
Figure GDA00035039715900001312
其中i＝Lk-1+1,Lk-1+2,…,Lk-1+Jk。新生粒子的分类可表示为：
Figure GDA00035039715900001313
步骤3.更新步骤：
(3.1)更新目标权值：
在时刻k接收到测量值之后，更新Lk-1+Jk个粒子权重:
Figure GDA0003503971590000141
其中，
Figure GDA0003503971590000142
其中，
Figure GDA0003503971590000143
表示排列系数，|Z|表示量测集的势，pK,k是杂波随机集的势分布，Zk\{z}表示在Zk中除了z的剩余量测。
Ξk(D,Z)＝{<D,Ψk,z(x)>:z∈Z}
Figure GDA0003503971590000144
其中，ej是一个j阶的初等对称函数。
ej({ρ1,ρ2,…,ρm})＝(-1)jαm-j/αm
其中，{ρ1,ρ2,…,ρm}是多项式αmxm+αm-1xm-1+...+α1x+α0的不同的根。
(3.2)更新势分布pk(n)：
Figure GDA0003503971590000145
(3.3)估计目标强度函数Dk：
Figure GDA0003503971590000146
其中，
Figure GDA0003503971590000147
是Parzen-Rosenblatt核函数，σd是核宽。
(3.3)量测更新后粒子的标记可表示为：
Figure GDA0003503971590000148
其中，i＝1,2,…Lk-1+Jk。更新后粒子的分类可以表示为：
Figure GDA0003503971590000149
步骤4.
计算当前时刻总的目标数目nk：
Figure GDA00035039715900001410
步骤5.更新时变参数粒子：
通过聚类获得同一目标簇的所有粒子，则时变参数粒子
Figure GDA0003503971590000151
的均值
Figure GDA0003503971590000152
和协方差Vk-1的估计可表示为：
Figure GDA0003503971590000153
其中，
Figure GDA0003503971590000154
是缩放因子，根据
Figure GDA0003503971590000155
的均值
Figure GDA0003503971590000156
和协方差Vk-1的估计，可以更新时变参数粒子
Figure GDA0003503971590000157
得到
Figure GDA0003503971590000158
步骤6.重采样步骤：
(6.1)通过权重对粒子集
Figure GDA0003503971590000159
进行重采样得到新的粒子集
Figure GDA00035039715900001510
重采样之后每个粒子被赋予相同的权重
Figure GDA00035039715900001511
其中，Lk＝Lk-1+Jk。
(6.2)重采样后的粒子与它们的父辈粒子具有相同的标记，如果：
Figure GDA00035039715900001512
分配标记：
Figure GDA00035039715900001513
步骤7.提取目标状态
(7.1)通过粒子聚类获得目标状态，并且聚类中心作为目标的估计状态：
Figure GDA00035039715900001514
其中，
Figure GDA00035039715900001515
Figure GDA00035039715900001516
是估计的目标数，round(·)表示舍入运算符。
(7.2)给新的目标分配新的标记
Figure GDA00035039715900001517
得到估计后的分类为
Figure GDA00035039715900001518
步骤8.航迹关联：
在前一时刻，我们得到的粒子分类集合为
Figure GDA00035039715900001519
在当前时刻，粒子的分类集合为
Figure GDA00035039715900001520
定义两个矩阵：
Figure GDA00035039715900001521
Figure GDA00035039715900001522
其中*表示满足条件的粒子数目，矩阵M表示k时刻，每一个分类中的粒子有多少是与k-1时刻的分类相对应。矩阵F说明在k时刻，每一个分类重采样的粒子有多少是与k-1时刻相对应。设定一个门限值ε1，如果对于前一时刻的目标
Figure GDA00035039715900001523
N为每个目标的粒子数，认为目标g消失。新生目标周围的粒子是从已经存在的目标粒子云中采样得到。在这种情况下，前一时刻的目标g在k时刻可能会分为两类。在矩阵M中这两个目标的粒子数目可能一样，所以需要看矩阵F。矩阵F指出有多少粒子被重采样，因为如果存活的目标跟踪准确时，从存活目标中重采样的粒子数要比从新生目标中重采样的粒子数要多。定义一个有效矩阵A，定义如下：
mg,h≥ε1N，Ag,h＝1
mg,h<ε1N，Ag,h＝0
航迹关联估计如下：
∑hAg,h＝0，删除目标表示Tk,g
∑hAg,h＝1，则目标g与h关联
∑hAg,h>1，取h＝argmaxFg,h，将目标g与h关联
步骤9.重复步骤2，继续跟踪下一时刻的多目标。
本发明的效果可通过以下实验进一步说明：
1.实验条件及参数
实验是在处理器为Intel Core i5-3470、3.2GHz，内存为4GB的Dell计算机平台上，采用MATLAB 2016仿真软件完成。模拟了二位跟踪场景，在二维跟踪场景下，测量数据在四个传感器中获得，四个传感器的位置分别是(0,0)m，(0,1×104)m，(1×104,0)m，以及(1×104，1×104)m。在时间k时刻，每个传感器输出所测量的接收信号的方位，由下式确定：
Figure GDA0003503971590000161
其中，
Figure GDA0003503971590000166
表示第i只传感器的位置，i＝1,2,3,4。wk是方差为
Figure GDA0003503971590000162
的零均值高斯白噪声。实验场景中有三个机动目标，目标1和目标2在整个模拟过程中保持活动状态，初始位置为(-3×103,5×103)m和(1.4×104,8×103)m。目标3在第10分钟产生，初始位置为(2×103,10.5×103)m，在第50分钟消失。三个目标的真实轨迹如图2所示，匀速运动(CV)模型和两个转弯(CT)模型的状态方程如下：
Figure GDA0003503971590000163
Figure GDA0003503971590000164
其中，
Figure GDA0003503971590000165
是第i个目标的状态向量。F(ω)是转弯模型的状态转移矩阵：
Figure GDA0003503971590000171
将两个转弯模型的转弯率设置为±9°/min，
Figure GDA0003503971590000172
是协方差为：
Figure GDA0003503971590000173
Figure GDA0003503971590000174
的零均值高斯白噪声。杂波被建模为观测空间上杂波率r＝10的泊松分布。目标的存活概率和被检测概率分别为PS＝0.99和Pd＝0.98。逆伽马分布的初始参数设置为a＝b＝1。假设目标采样粒子的最大数目为1500，最小数目为300。采用目标数目估计和最佳子模式分配(OSPA)距离作为提出方法评价的性能指标，其中，OSPA距离定义为
Figure GDA0003503971590000175
其中，X＝{x1,x2,…xm}，Y＝{y1,y2,…yn}是任意有限的子集，1≤P<∞,c>0,m,n∈N0＝{0,1,2,…}。如果
Figure GDA0003503971590000176
OSPA距离参数P＝2，c＝1000。实验结果为100次蒙特卡罗仿真的统计结果。
2.实验及结果分析
本发明方法的实验，主要从以下三个方面开展：
实验1：PF-APE-CPHD与PF-APE-PHD，PF-MM-PHD，PF-MM-CPHD算法性能对比
本组实验场景中，对于PF-APE-CPHD，PF-APE-PHD，PF-MM-PHD，PF-MM-CPHD,已知测量噪声的标准偏差设置为σ＝0.03，对于PF-APE-CPHD，PF-APE-PHD，PF-MM-PHD作为时变参数的角速度ω是未知的，以验证提出算法对时变参数的自适应能力。对于PF-MM-CPHD,采用转弯率ω＝9°/min和ω＝-9°/min的转弯模型。对于PF-APE-PHD，PF-MM-PHD，PF-MM-CPHD新生目标强度是已知的，采用泊松分布对新生目标的出生过程进行建模，即：
Figure GDA0003503971590000181
其中，i＝1,2,3。
Figure GDA0003503971590000182
Figure GDA0003503971590000183
Figure GDA0003503971590000184
Figure GDA0003503971590000185
以验证提出算法对新生目标自动识别的能力。
图3给出了本发明方法及PF-APE-PHD，PF-MM-PHD，PF-MM-CPHD的跟踪OSPA距离，很明显，本发明方法优于PF-MM-PHD及PF-MM-CPHD算法，和PF-APE-PHD算法性能接近。这也表明所提出的方法可以适应目标机动参数的变化，并且在新生目标强度未知的情况下，自动识别新生目标，达到和PF-APE-PHD算法性能接近的水平。
图4显示了通过PF-APE-CPHD，PF-APE-PHD，PF-MM-PHD和PF-MM-CPHD过滤器获得的平均目标数估计。可以看出，所提出的PF-APE-CPHD算法明显优于PF-MM-PHD和PF-MM-CPHD的目标数量估计，和PF-APE-PHD算法目标数量估计接近。因为提出算法可以有效地联合估计未知模型参数ω，可以较好地调整目标的运动模型，而对于PF-MM-PHD,PF-MM-CPHD算法，由于多个模型的相互作用，模型干扰会影响精度，这是基于IMM技术的不可避免的现象。
实验2：具有未知测量噪声方差的多个突然机动目标跟踪实验
本组实验场景中，测量噪声的真实标准偏差固定在σ＝0.01rad，但是对于本发明方法是未知的，PF-APE-CPHD方法将其与时变转弯率ω和目标状态进行联合估计。为了便于分析提出方法的性能，同时假设在PF-APE-CPHD滤波器中采用不同的已知量测噪声方差进行滤波，量测噪声方差分别取σ＝0.005,0.01,0.03,0.05,0.06。
图5和图6给出了本发明方法在具有未知测量噪声方差情况下多个机动目标跟踪的平均OSPA距离和目标数估计。可以看出，当测量噪声方差与目标状态联合估计时，PF-APE-CPHD算法的性能接近于采用已知真实量测噪声方差σ＝0.01rad的估计性能，因此，提出方法能够近似估计未知的噪声参数，适应对未知噪声参数的复杂场景下目标跟踪。此外，从图中还可以看出，以不正确的测量噪声方差进行滤波时，算法跟踪结果会显著下降，也即不正确的噪声参数会导致目标的运动模型不匹配，影响跟踪性能。
实验3：机动多目标航迹维持
本组实验场景中，对于PF-APE-CPHD，已知测量噪声的标准偏差设置为σ＝0.01，对于PF-APE-CPHD，时变参数(角速度ω)未知，验证提出方法对航迹跟踪性能。
图7给出了在PF-APE-CPHD算法中，运用粒子标识技术可以有效跟踪每个目标的航迹，由于机动目标的跟踪存在一定的误差和受到粒子标识中阈值的影响，在某些时刻的的航迹跟踪存在一定的偏差。

Claims (9)

1.一种未知新生强度的多目标跟踪方法，其特征在于，所述方法是将采样的目标粒子按照属于各自目标进行分类，获得粒子集
Figure FDA0003503971580000011
其中，k表示当前时刻，
Figure FDA0003503971580000012
表示分类数目；选择
Figure FDA0003503971580000013
中同类的所有粒子，计算每一类粒子的位置均值
Figure FDA0003503971580000014
计算每一类粒子位置均值与当前时刻每个量测位置的距离
Figure FDA0003503971580000015
其中
Figure FDA0003503971580000016
Mk表示当前时刻量测的个数；根据距离最短的
Figure FDA0003503971580000017
进行存活目标判断，表明第i类目标粒子与当前的第j个量测对应，k时刻未配对的量测记为Λk；计算每一类粒子位置均值与k+1时刻每个量测位置的距离
Figure FDA0003503971580000018
Mk+1表示k+1时刻量测的个数，根据距离最短的
Figure FDA0003503971580000019
进行判断存活目标判断，如果距离
Figure FDA00035039715800000110
小于阈值U，表明第i类目标粒子与当前的第j个量测对应，k+1时刻未配对的量测记为Λk+1，对于k时刻剩下的量测，计算与k+1时刻剩下量测位置的距离
Figure FDA00035039715800000111
其中，
Figure FDA00035039715800000112
表示k时刻量测集中没有和目标匹配的量测个数，
Figure FDA00035039715800000113
表示k+1时刻量测集中没有和目标匹配的量测的个数，若距离
Figure FDA00035039715800000114
小于阈值U，则k时刻第i个未关联的量测为新生目标，并根据该量测采样N个新生目标粒子，k时刻如果没有被关联的量测，判断为杂波；
所述的阈值U是指根据多个目标的最大运动速度和采样间隔的乘积设定最大距离反映当前时刻量测与下一时刻量测关联情况下的最大距离，当小于最大距离代表有新生目标生成。
2.一种机动多目标跟踪或航迹管理方法，所述方法包括如下步骤：
(1)初始化步骤：
(1a)初始时刻，假设有n0个目标，每个目标采样N个粒子，总粒子数为L0＝N×n0；初始目标状态集为X0，目标存在概率为Ps，目标检测概率为Pd，势分布为p0；
(1b)初始化采样粒子并对粒子进行标识，设P0(X0,θ0,ξ0)为先验的联合概率分布；从P0(X0,θ0,ξ0)中采用初始粒子集
Figure FDA00035039715800000115
每个粒子的权重设为
Figure FDA00035039715800000116
初始时刻，粒子标记表示为
Figure FDA00035039715800000117
初始时刻粒子的分类可表示为
Figure FDA00035039715800000118
Figure FDA00035039715800000119
表示分类数目；
(2)当k≥1时，产生2Lk-1个粒子，其中，从后验分布
Figure FDA00035039715800000120
中产生Lk-1个预测粒子
Figure FDA00035039715800000121
权重为
Figure FDA00035039715800000122
其中，i＝1,2,…,Lk-1，
Figure FDA00035039715800000123
是时变参数粒子
Figure FDA00035039715800000124
的均值，Zk是量测集；然后从从后验分布
Figure FDA00035039715800000125
中产生其余的Lk-1个粒子，权重为
Figure FDA0003503971580000021
其中，
Figure FDA0003503971580000022
为是从时变参数服从的先验分布Pξ(ξ0)中采样获得，i＝Lk-1+1,Lk-1+2,…,2Lk-1；
(3)计算量测噪声协方差：
(3a)假设
Figure FDA0003503971580000023
为静态参数粒子
Figure FDA0003503971580000024
的充分统计量，参数粒子通过
Figure FDA0003503971580000025
获得，其中～表示采样，P()表示分布，
Figure FDA0003503971580000026
表示粒子
Figure FDA0003503971580000027
从分布
Figure FDA0003503971580000028
中采样获得；考虑静态未知参数Vk为未知的量测噪声方差，其共轭先验可采用参数为a和b的逆伽马分布IG(a,b)近似；设定逆伽马分布参数，
Figure FDA0003503971580000029
Figure FDA00035039715800000210
其中，l＝1,…,d，d为量测噪声协方差R的维度；
(3b)计算量测噪声协方差
Figure FDA00035039715800000211
其中，n＝1,…,Lk-1为迭代次序，Lk-1为最大迭代次数；
(3c)如果n≤N，更新参数
Figure FDA00035039715800000212
返回步骤(3b)；
(4)计算每个粒子量测噪声协方差
Figure FDA00035039715800000213
利用量测噪声协方差计算前Lk-1个粒子对应于时变参数无变化点
Figure FDA00035039715800000214
预测似然成比例的另一个权重
Figure FDA00035039715800000215
其中i＝1,2,,…Lk-1；计算剩余Lk-1个粒子对应于时变参数变化点
Figure FDA00035039715800000216
的预测似然成比例的另一个权重
Figure FDA00035039715800000217
其中i＝Lk-1+1,Lk-1+2,…,2Lk-1；
(5)选择
Figure FDA00035039715800000218
中同类的所有粒子，计算每一类粒子的位置均值
Figure FDA00035039715800000219
计算每一类粒子位置均值与当前时刻每个量测位置的距离
Figure FDA00035039715800000220
其中
Figure FDA00035039715800000221
Mk表示当前时刻量测的个数；根据距离最短的
Figure FDA00035039715800000222
进行判断存活目标判断，表明第i类目标粒子与当前的第j个量测对应，k时刻未配对的量测记为Λk；计算每一类粒子位置均值与k+1时刻每个量测位置的距离
Figure FDA00035039715800000223
Mk+1表示k+1时刻量测的个数，根据距离最短的
Figure FDA00035039715800000224
进行判断存活目标判断，如果距离
Figure FDA00035039715800000225
小于阈值U，表明第i类目标粒子与当前的第j个量测对应，k+1时刻未配对的量测记为Λk+1，对于k时刻剩下的量测，计算与k+1时刻剩下量测位置的距离
Figure FDA00035039715800000226
其中，
Figure FDA00035039715800000227
表示k时刻量测集中没有和目标匹配的量测个数，
Figure FDA00035039715800000228
表示k+1时刻量测集中没有和目标匹配的量测的个数，若距离
Figure FDA00035039715800000229
小于阈值U，则k时刻第i个未关联的量测为新生目标，则根据当前对应量测位置采样N个新生目标粒子；
(6)选择2Lk-1个粒子中的Lk-1个粒子，将它们的索引表示为：l(i)∈{1，2，...，2Lk-1}，其中，i＝1，2，...，Lk-1；对于i＝1，2，...，Lk-1，选择索引li来自{1，2，...，Lk-1}的概率为
Figure FDA0003503971580000031
来自{Lk-1+1，Lk-1+2，...，2Lk-1}的概率为
Figure FDA0003503971580000032
其中β是发生突变的概率，并假设已知；
(7)预测势分布pk|k-1；
(8)预测粒子标记
Figure FDA0003503971580000033
和粒子分类
Figure FDA0003503971580000034
(9)更新目标权重
Figure FDA0003503971580000035
得到
Figure FDA0003503971580000036
(10)更新势分布pk；
(11)估计目标强度函数Dk；
(12)更新粒子标记
Figure FDA0003503971580000037
和粒子分类
Figure FDA0003503971580000038
(13)计算当前时刻总的目标数目nk；
(14)聚类获得同一目标簇的粒子，计算时变参数粒子
Figure FDA0003503971580000039
的均值
Figure FDA00035039715800000310
和协方差Vk-1的估计，更新时变参数
Figure FDA00035039715800000311
的值，得到
Figure FDA00035039715800000312
(15)重采样：
(15a)通过权重对粒子集
Figure FDA00035039715800000313
进行重采样得到新的粒子集
Figure FDA00035039715800000314
重采样之后每个粒子被赋予相同的权重
Figure FDA00035039715800000315
其中，Lk＝Lk-1+Jk，Jk是k时刻新生目标粒子数；
(15b)重采样后的粒子与它们的父辈粒子具有相同的标记，对重采样后的粒子进行标记，得到
Figure FDA00035039715800000316
其中，j＝1，2，...，Lk-1+Jk；
(16)通过聚类粒子获得目标状态，并且聚类中心是估计的状态
Figure FDA00035039715800000317
给新的分类分配新的标记
Figure FDA00035039715800000318
得到估计后的分类为
Figure FDA00035039715800000319
(17)通过粒子标记进行航迹关联得到每个目标的航迹，若下一时刻观测信息到达，转到步骤(2)进行迭代；否则，目标跟踪过程结束。
3.根据权利要求2所述的方法，其特征在于，所述步骤(4)中利用量测噪声方差计算前Lk-1个粒子对应于时变参数无变化点
Figure FDA00035039715800000320
预测似然成比例的另一个权重
Figure FDA00035039715800000321
Figure FDA00035039715800000322
其中，i＝1，2，...，Lk-1，计算剩余Lk-1个粒子对应于时变参数变化点
Figure FDA00035039715800000323
的预测似然成比例的另一个权重
Figure FDA00035039715800000324
Figure FDA0003503971580000041
其中，i＝Lk-1+1，Lk-1+2，...，2Lk-1。
4.根据权利要求2所述的方法，其特征在于，所述步骤(7)中计算预测势分布pk|k-1的公式为：
Figure FDA0003503971580000042
Figure FDA0003503971580000043
其中，pΓ，k是新生目标随机集的势分布，<.，.>表示内积运算；
Figure FDA00035039715800000410
表示二项式系数。
5.根据权利要求2所述的方法，其特征在于，所述步骤(9)中更新目标权重
Figure FDA0003503971580000044
得到
Figure FDA0003503971580000045
涉及的公式为：
在时刻k接收到测量值之后，更新Lk-1+Jk个粒子权重：
Figure FDA0003503971580000046
其中，
Figure FDA0003503971580000047
其中，
Figure FDA0003503971580000048
表示排列系数，|Z|表示量测集的势，pK，k是杂波随机集的势分布，Zk\{z}表示在Zk中除了z的剩余量测；
Ξk(D，Z)＝{<D，Ψk，z(x)>：z∈Z}
Figure FDA0003503971580000049
其中，ej是一个j阶的初等对称函数：
ej({ρ1，ρ2，...，ρm})＝(-1)jαm-j/αm
其中，{ρ1，ρ2，...，ρm}是多项式αmxm+αm-1xm-1+...+α1x+α0的不同的根。
6.根据权利要求2所述的方法，其特征在于，所述步骤(10)中更新势分布pk(n)的计算公式为：
Figure FDA0003503971580000051
7.根据权利要求2所述的方法，其特征在于，所述步骤(11)中估计目标强度函数Dk：
Figure FDA0003503971580000052
其中，
Figure FDA0003503971580000053
是ParZen-Rosenblatt核函数，σd是核宽。
8.根据权利要求2所述的方法，其特征在于，所述步骤(12)中量测更新后粒子的标记可表示为：
Figure FDA0003503971580000054
其中，i＝1，2，...Lk-1+Jk。
9.根据权利要求2所述的方法，其特征在于，所述步骤(12)中更新后粒子的分类可以表示为：
Figure FDA0003503971580000055
CN202010466368.2A 2020-05-28 2020-05-28 一种未知新生强度的机动多目标跟踪与航迹维持方法 Active CN111562571B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010466368.2A CN111562571B (zh) 2020-05-28 2020-05-28 一种未知新生强度的机动多目标跟踪与航迹维持方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010466368.2A CN111562571B (zh) 2020-05-28 2020-05-28 一种未知新生强度的机动多目标跟踪与航迹维持方法

Publications (2)

Publication Number Publication Date
CN111562571A CN111562571A (zh) 2020-08-21
CN111562571B true CN111562571B (zh) 2022-04-29

Family

ID=72069795

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010466368.2A Active CN111562571B (zh) 2020-05-28 2020-05-28 一种未知新生强度的机动多目标跟踪与航迹维持方法

Country Status (1)

Country Link
CN (1) CN111562571B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111562571B (zh) * 2020-05-28 2022-04-29 江南大学 一种未知新生强度的机动多目标跟踪与航迹维持方法
CN112215146B (zh) * 2020-10-12 2022-04-22 西安交通大学 基于随机有限集的弱小目标联合检测与跟踪系统及方法
CN113987980B (zh) * 2021-09-23 2022-05-20 北京连山科技股份有限公司 一种图形物理phd流行仿真实现方法
CN114325686B (zh) * 2021-12-23 2024-05-03 中国人民解放军国防科技大学 基于smc-phd滤波器的多目标跟踪方法
CN115015906A (zh) * 2022-05-31 2022-09-06 西安电子科技大学 一种目标航迹估计方法以及相关装置

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103345577A (zh) * 2013-06-27 2013-10-09 江南大学 变分贝叶斯概率假设密度多目标跟踪方法
CN105761276A (zh) * 2015-12-15 2016-07-13 江南大学 基于迭代ransac自适应新生目标强度估计的gm-phd多目标跟踪算法
CN106526585A (zh) * 2016-10-26 2017-03-22 中国人民解放军空军工程大学 基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法
CN106772251A (zh) * 2016-11-18 2017-05-31 中国船舶重工集团公司第七二四研究所 一种时差定位系统多波束资源优化调度方法
CN107037423A (zh) * 2016-11-09 2017-08-11 谭顺成 结合幅值信息的phd滤波多目标跟踪方法
CN110376582A (zh) * 2019-01-24 2019-10-25 西安电子科技大学 自适应gm-phd的机动目标跟踪方法
CN110619409A (zh) * 2018-06-19 2019-12-27 新智数字科技有限公司 一种自适应扰动量子粒子群的泛能站调度方法及装置
CN110780269A (zh) * 2019-10-08 2020-02-11 河海大学 自适应新生强度下基于gm-phd滤波器的显式多目标跟踪方法
CN111562571A (zh) * 2020-05-28 2020-08-21 江南大学 一种未知新生强度的机动多目标跟踪与航迹维持方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5862023B2 (ja) * 2011-03-04 2016-02-16 日本電気株式会社 目標追跡システム及び目標追跡方法
US9562965B2 (en) * 2011-05-04 2017-02-07 Jacques Georgy Two-stage filtering based method for multiple target tracking

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103345577A (zh) * 2013-06-27 2013-10-09 江南大学 变分贝叶斯概率假设密度多目标跟踪方法
CN105761276A (zh) * 2015-12-15 2016-07-13 江南大学 基于迭代ransac自适应新生目标强度估计的gm-phd多目标跟踪算法
CN106526585A (zh) * 2016-10-26 2017-03-22 中国人民解放军空军工程大学 基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法
CN107037423A (zh) * 2016-11-09 2017-08-11 谭顺成 结合幅值信息的phd滤波多目标跟踪方法
CN106772251A (zh) * 2016-11-18 2017-05-31 中国船舶重工集团公司第七二四研究所 一种时差定位系统多波束资源优化调度方法
CN110619409A (zh) * 2018-06-19 2019-12-27 新智数字科技有限公司 一种自适应扰动量子粒子群的泛能站调度方法及装置
CN110376582A (zh) * 2019-01-24 2019-10-25 西安电子科技大学 自适应gm-phd的机动目标跟踪方法
CN110780269A (zh) * 2019-10-08 2020-02-11 河海大学 自适应新生强度下基于gm-phd滤波器的显式多目标跟踪方法
CN111562571A (zh) * 2020-05-28 2020-08-21 江南大学 一种未知新生强度的机动多目标跟踪与航迹维持方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"Multiple-model probability hypothesis density filter for tracking maneuvering targets";Punithakumar K 等;《 IEEE Aerospace and Electronic Systems Magazine》;20081231;第87-98页 *
"Multi-sensor GIW-PHD filter for multiple extended target tracking";Peng Li 等;《The 27th Chinese Control and Decision Conference (2015 CCDC)》;20151231;第5620-5625页 *

Also Published As

Publication number Publication date
CN111562571A (zh) 2020-08-21

Similar Documents

Publication Publication Date Title
CN111562571B (zh) 一种未知新生强度的机动多目标跟踪与航迹维持方法
CN110084831B (zh) 基于YOLOv3多伯努利视频多目标检测跟踪方法
CN110503071B (zh) 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法
CN106572493B (zh) Lte网络中的异常值检测方法及系统
CN110084836B (zh) 基于深度卷积特征分层响应融合的目标跟踪方法
CN109508444B (zh) 区间量测下交互式多模广义标签多伯努利的快速跟踪方法
CN107831490A (zh) 一种改进的多扩展目标跟踪方法
CN110501671B (zh) 一种基于测量分配的目标跟踪方法及装置
CN105405151A (zh) 基于粒子滤波和加权Surf的抗遮挡目标跟踪方法
CN105354860B (zh) 基于箱粒子滤波的扩展目标CBMeMBer跟踪方法
CN105182291A (zh) 自适应目标新生强度的phd平滑器的多目标跟踪方法
CN107871156B (zh) 基于信息素预测的蚁群多细胞跟踪系统
CN109214432B (zh) 一种多传感器多目标联合检测、跟踪与分类方法
CN111830501B (zh) Hrrp历史特征辅助的信号模糊数据关联方法及系统
CN111488552B (zh) 基于高斯混合概率假设密度的紧邻多目标跟踪方法
CN110889862A (zh) 一种网络传输攻击环境中多目标跟踪的组合测量方法
CN114137526A (zh) 基于标签的车载毫米波雷达多目标检测方法和系统
CN115204212A (zh) 一种基于stm-pmbm滤波算法的多目标跟踪方法
CN113406623A (zh) 基于雷达高分辨距离像的目标识别方法、装置及介质
CN111711432A (zh) 一种基于ukf和pf混合滤波的目标跟踪算法
CN109671096B (zh) 一种时空近邻目标检测及网格聚类量测划分下的多扩展目标跟踪方法
CN113409363B (zh) 一种基于bp-pmbm滤波算法的多目标跟踪方法
CN113537411A (zh) 一种基于毫米波雷达的改进模糊聚类方法
CN111274529B (zh) 一种鲁棒的高斯逆威沙特phd多扩展目标跟踪算法
CN115619825A (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
GR01 Patent grant
GR01 Patent grant