CN112906743B - 一种快速多传感器集势概率假设密度滤波方法 - Google Patents

一种快速多传感器集势概率假设密度滤波方法 Download PDF

Info

Publication number
CN112906743B
CN112906743B CN202110069747.2A CN202110069747A CN112906743B CN 112906743 B CN112906743 B CN 112906743B CN 202110069747 A CN202110069747 A CN 202110069747A CN 112906743 B CN112906743 B CN 112906743B
Authority
CN
China
Prior art keywords
measurement
sensor
hypothesis density
probability hypothesis
matrix
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
CN202110069747.2A
Other languages
English (en)
Other versions
CN112906743A (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202110069747.2A priority Critical patent/CN112906743B/zh
Publication of CN112906743A publication Critical patent/CN112906743A/zh
Application granted granted Critical
Publication of CN112906743B publication Critical patent/CN112906743B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/251—Fusion techniques of input or preprocessed data
    • 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
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Algebra (AREA)
  • Probability & Statistics with Applications (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明属于雷达跟踪领域，目的在于提供一种快速多传感器集势概率假设密度滤波方法，并使用高斯混合技术进行实现。该方法通过构建简化处理后的权重矩阵和代价矩阵，减少了后续量测分组中可能的分组数，避免了分组过程中可能出现的重复问题，极大地提高了算法效率，同时也提升了跟踪性能，获得较传统方法更加稳定准确的目标状态和目标数估计，具有很好的工程应用前景。

Description

一种快速多传感器集势概率假设密度滤波方法
技术领域
本发明涉及目标跟踪技术领域，具体地，涉及一种快速多传感器集势概率假设密度(Fast Multisensor Cardinalized Probability Hypothesis Density，FMS-CPHD)滤波方法。
背景技术
随着目标跟踪场景越发复杂，对目标跟踪技术的要求也在不断提升，单传感器性能往往无法满足需求，多传感器协同跟踪成为当前监视场景的普遍需求。相对于单传感器跟踪，多传感器跟踪最大的问题就是耗时长，时效性差，当跟踪环境比较复杂时可能无法满足实时性要求。同时，多传感器跟踪还存在跟踪性能不稳定，跟踪性能不能保证比其中某一个传感器性能更好等情况。
集中式多传感器融合跟踪可以保证跟踪性能较任意单传感器更好，跟踪性能有保证，但是存在计算复杂度高，计算效率低的情况。基于随机有限集(Random Finite Set，RFS)理论的多传感器集势概率假设密度(Multisensor Cardinalized ProbabilityHypothesis Density，MS-CPHD)滤波器是一种集中式多传感器融合跟踪方法，具有较好的跟踪性能，但是计算复杂度太高，实时性差，而且存在性能不稳定的情况，因此很难在工程实践中获得应用。
发明内容
本发明的目的在于提供一种快速多传感器集势概率假设密度(FMS-CPHD)滤波方法，并使用高斯混合技术进行实现。该方法通过构建简化处理后的权重矩阵和代价矩阵，减少了后续量测分组中可能的分组数，避免了分组过程中出现的重复问题，极大地提高了算法效率，同时也提升了跟踪性能，获得较传统方法更加稳定准确的目标状态和目标数估计。
实现本发明目的的技术解决方案为：一种快速多传感器集势概率假设密度滤波方法，包括以下步骤：
先给出统一的符号定义：(·)k表示k时刻的值，则(·)0表示初始值，(·)k|k1表示k-1时刻对k时刻的预测值，(·)(i)表示与索引号为i的高斯分量有关的物理量，(·)γ表示与传感器γ有关的物理量，(·)l表示与目标航迹l有关的物理量；
再给出高斯分量的定义：
Figure BDA0002905571730000011
表示k时刻索引号为i的高斯分量，其权重为
Figure BDA0002905571730000012
均值为
Figure BDA0002905571730000013
协方差为
Figure BDA0002905571730000014
被标签
Figure BDA0002905571730000015
标记；索引号i∈{1，2，…，Ik}，其中Ik表示k时刻高斯分量的数量；标签
Figure BDA0002905571730000016
其中
Figure BDA0002905571730000017
表示k时刻的标签数(目标数)，被同一个标签标记的高斯分量认为来自同一个目标航迹；
S1，场景中总共有Υ个传感器，传感器编号γ∈{1，2，...，}；设初始时刻k＝0，初始化势分布(目标数分布)ρ0(n)和概率假设密度D0(x)，具体为：
初始化势分布为ρ0(n)，其中n表示目标数为n且n∈{0，1，…，Nmax}，Nmax为最大可能目标数；ρ0(n)可以根据情况选择合适分布，这里选择使用二项分布：
Figure BDA0002905571730000021
其中C表示组合数符号，
Figure BDA0002905571730000022
r为目标出现概率；
初始化概率假设密度D0(x)的形式为：
Figure BDA0002905571730000023
S2，已有上一时刻(k-1时刻)的概率假设密度Dk-1(x)和势分布ρk-1(n)，对当前时刻(k时刻)进行预测，得到预测概率假设密度Dk|k-1(x)和预测势分布ρk|k-1(n)；包括以下步骤：
上一时刻的概率假设密度Dk-1(x)的形式为：
Figure BDA0002905571730000024
对当前时刻的势分布进行预测为：
Figure BDA0002905571730000025
其中ρB(·)是新目标的势分布，pS是目标存活概率，s和t表示两个整数；
对当前时刻的概率假设密度进行预测为：
Dk|k-1(x)＝DS，k|k-1(x)+DB，k|k-1(x)
其中DS，k|k-1(x)表示存活目标的概率假设密度，DB，k|k-1(x)表示新目标概率假设密度；
DS，k|k-1(x)由下式给出：
Figure BDA0002905571730000026
其中
Figure BDA0002905571730000027
是存活目标标签且
Figure BDA0002905571730000028
和
Figure BDA0002905571730000029
分别是存活目标高斯分量的均值和协方差，由下两式分别计算：
Figure BDA00029055717300000210
Figure BDA00029055717300000211
其中F和Q分别表示状态转移矩阵和过程噪声协方差矩阵，T表示矩阵转置；
DB，k|k-1(x)由下式给出：
Figure BDA0002905571730000031
其中
Figure BDA0002905571730000032
和
Figure BDA0002905571730000033
分别是新目标高斯分量的权重，均值和协方差，IB，k|k-1是新目标高斯分量数，新目标标签
Figure BDA0002905571730000034
是新目标数；
则预测概率假设密度Dk|k-1(x)为：
Figure BDA0002905571730000035
其中Ik|k-1＝Ik-1+IB，k|k-1，标签
Figure BDA0002905571730000036
S3，获得传感器γ当前时刻的量测集合
Figure BDA0002905571730000037
其中
Figure BDA0002905571730000038
是量测状态，
Figure BDA0002905571730000039
是量测数量，通过各传感器的量测集合，计算各传感器的权重矩阵和代价矩阵；包括以下步骤：
S3.1传感器γ的代价矩阵用
Figure BDA00029055717300000310
表示，矩阵第l行，第j列的元素
Figure BDA00029055717300000311
表示航迹l与量测
Figure BDA00029055717300000312
的关联代价，其中
Figure BDA00029055717300000313
由下式计算：
Figure BDA00029055717300000314
其中：
Figure BDA00029055717300000315
表示标记为l的高斯分量索引号的集合；
Figure BDA00029055717300000316
H是量测矩阵，
Figure BDA00029055717300000317
是量测噪声协方差矩阵；
Figure BDA00029055717300000318
表示传感器γ的检测概率；
κ(·)表示杂波的强度函数；
S3.2传感器γ的权重矩阵由
Figure BDA00029055717300000319
表示，其大小与
Figure BDA00029055717300000320
相同，第l行，第j列的元素
Figure BDA00029055717300000321
由下式计算：
Figure BDA00029055717300000322
S3.3计算量测
Figure BDA00029055717300000323
的更新权重：
Figure BDA0002905571730000041
设置量测剔除阈值θz，量测剔除阈值一般可设置为不大于10-3的数，通常取θz＝10-3，如果
Figure BDA0002905571730000042
则认为
Figure BDA0002905571730000043
是杂波，将其所在列从代价矩阵
Figure BDA0002905571730000044
中剔除；剔除后矩阵
Figure BDA0002905571730000045
剩余
Figure BDA0002905571730000046
列；
S3.4计算传感器γ代价矩阵的行列差：
Figure BDA0002905571730000047
从各传感器代价矩阵的行列差
Figure BDA0002905571730000048
中取最大值为
Figure BDA0002905571730000049
如果Δk＞0，则将Δk+1个代价矩阵
Figure BDA00029055717300000410
纵向拼接为新的代价矩阵：
Figure BDA00029055717300000411
S3.5代价矩阵
Figure BDA00029055717300000412
的关联矩阵用
Figure BDA00029055717300000413
表示，其大小与
Figure BDA00029055717300000414
相同，其中第l行，第j列的元素为al，j，al，j等于0或1，且每行每列最多只有1个al，j等于1；则关联矩阵
Figure BDA00029055717300000415
的代价为
Figure BDA00029055717300000416
其中
Figure BDA00029055717300000417
表示此时
Figure BDA00029055717300000418
的行数，
Figure BDA00029055717300000419
表示
Figure BDA00029055717300000420
中第l行，第j列的元素；
通过寻优算法寻找代价最小的关联矩阵
Figure BDA00029055717300000421
例如匈牙利寻优算法。然后求矩阵
Figure BDA00029055717300000422
与
Figure BDA00029055717300000423
的Hadamard积，得到最终的代价矩阵
Figure BDA00029055717300000424
其中符号*表示Hadamard积运算；
S4，利用S3获得的代价矩阵
Figure BDA00029055717300000425
通过航迹对各传感器的量测进行分组，然后用各航迹的量测分组组成全局量测划分，包括以下步骤：
S4.1用每个航迹对各传感器的量测进行分组，具体过程为：
对于航迹l，依次从传感器γ∈{1，2，...，Υ}的量测集合
Figure BDA00029055717300000426
中选择一个量测，且只有在传感器γ的代价矩阵
Figure BDA00029055717300000427
中
Figure BDA00029055717300000428
的量测可以被选择，不选择量测则表示传感器γ漏检该目标航迹l；完成后将所有选择的量测组成航迹l的第1个量测分组
Figure BDA00029055717300000429
S4.2重复S4.1的过程，找到航迹l所有可能的量测分组
Figure BDA00029055717300000430
其中Ll是航迹l的量测分组数；
S4.3对所有的航迹进行S4.1～S4.2的过程，找到所有航迹的量测分组；
S4.4然后用各航迹的量测分组组成全局量测划分，具体过程为：
依次从航迹
Figure BDA0002905571730000051
的量测分组中选择一个分组
Figure BDA0002905571730000052
其中
Figure BDA0002905571730000053
将这些量测分组组成一个全局量测划分
Figure BDA0002905571730000054
其中
Figure BDA0002905571730000055
表示来自所有传感器的杂波组成的集合；用
Figure BDA0002905571730000056
表示量测分组
Figure BDA0002905571730000057
中来自传感器γ的量测数量，则ψ中来自传感器γ的量测数量表示为
Figure BDA0002905571730000058
ψ中总的量测数量表示为
Figure BDA0002905571730000059
S4.5重复该过程，找到所有可能的量测划分，所有量测划分组成集合Ψ；
S5，通过S2获得的预测概率假设密度Dk|k-1(x)和预测势分布ρk|k-1(n)，以及S4获得的量测划分集合Ψ进行更新，获得更新的概率假设密度Dk(x)和更新势分布ρk(n)；包括以下步骤：
S5.1更新的势分布ρk(n)由下式计算：
Figure BDA00029055717300000510
其中：
Figure BDA00029055717300000511
Figure BDA00029055717300000512
是传感器γ的杂波的势分布；
Figure BDA00029055717300000513
Figure BDA00029055717300000514
Figure BDA00029055717300000515
P表示排列数符号，
Figure BDA00029055717300000516
Figure BDA00029055717300000517
表示量测分组
Figure BDA00029055717300000518
的似然比，由下式计算：
Figure BDA00029055717300000519
其中：
Figure BDA0002905571730000061
cγ(·)是杂波的空间分布函数；
S5.2索引号为i的高斯分量在量测分组
Figure BDA0002905571730000062
中的归一化伪似然
Figure BDA0002905571730000063
由下式计算：
Figure BDA0002905571730000064
S5.3更新的概率假设密度Dk(x)为：
Dk(x)＝DE，k(x)+DU，k(x)
其中DE，k(x)和DU，k(x)分别表示遗留概率假设密度和量测更新概率假设密度，分别为：
Figure BDA0002905571730000065
Figure BDA0002905571730000066
其中：
Figure BDA0002905571730000067
而α0由下式计算：
Figure BDA0002905571730000068
Figure BDA0002905571730000069
表示全局量测划分ψ中量测分组
Figure BDA00029055717300000610
更新的概率假设密度，由下式给出：
Figure BDA00029055717300000611
其中比例系数
Figure BDA00029055717300000612
由下式计算：
Figure BDA00029055717300000613
量测分组
Figure BDA0002905571730000071
更新的高斯分量权重
Figure BDA0002905571730000072
由下式计算：
Figure BDA0002905571730000073
其中αψ由下式计算：
Figure BDA0002905571730000074
Figure BDA0002905571730000075
是量测
Figure BDA0002905571730000076
更新的高斯分量均值，
Figure BDA0002905571730000077
是更新协方差，分别由下式计算：
Figure BDA0002905571730000078
Figure BDA0002905571730000079
Figure BDA00029055717300000710
则量测更新概率假设密度DU，k(x)为：
Figure BDA00029055717300000711
其中
Figure BDA00029055717300000712
新标签
Figure BDA00029055717300000713
由量测分组
Figure BDA00029055717300000714
决定，
Figure BDA00029055717300000715
是航迹l选择的量测分组，则
Figure BDA00029055717300000716
(5.4)重新整理DU，k(x)中高斯分量的索引号，可将DU，k(x)整理为：
Figure BDA00029055717300000717
其中
Figure BDA00029055717300000718
为所有ψ中量测数量之和；
将DE，k(x)重新整理为：
Figure BDA00029055717300000719
其中
Figure BDA00029055717300000720
最终，所有高斯分量组成当前时刻更新的概率假设密度Dk(x)：
Figure BDA00029055717300000721
其中Ik＝Ik|k-1+Ik|k-1|Ψ|；
S6，对高斯分量进行剪枝与合并，估计目标数和目标状态，包括以下步骤：
S6.1设置高斯分量剔除阈值wP，高斯分量剔除阈值wP一般可设置为不大于10-5的数，通常取wP＝10-5，如果
Figure BDA0002905571730000081
将权重
Figure BDA0002905571730000082
对应的高斯分量删除；
S6.2设置合并距离阈值dM，合并距离阈值dM一般根据跟踪场景的大小决定，通常取dM＝4m，如果具有相同标签的高斯分量之间的距离小于dM，则将这些高斯分量合并；
S6.3设置高斯分量总数阈值Imax，高斯分量总数阈值Imax一般可设置为不小于100的整数，通常取Imax＝200，如果高斯分量数目大于200，则将权重排在200之后的高斯分量剔除；
S6.4目标数估计
Figure BDA0002905571730000083
就是使更新势分布取最大值的整数，即：
Figure BDA0002905571730000084
S6.5航迹l的高斯分量权重和为：
Figure BDA0002905571730000085
设置航迹剔除阈值θl，航迹剔除阈值θl一般可设置为不大于10-3的数，通常取θl＝10-4，如果航迹l的高斯分量权重w(l)θl，则将索引号
Figure BDA0002905571730000086
的高斯分量剔除；
将航迹高斯分量权重和按从大到小排序，取前
Figure BDA0002905571730000087
条航迹，然后从每个航迹对应的高斯分中选择权重最大的高斯分量，其均值就是目标状态的估计，比如航迹l的状态估计为：
Figure BDA0002905571730000088
S7，重复S2～S6，直到不再需要跟踪为止。
本发明具有以下技术效果：
1.本发明提出的FMS-CPHD滤波方法在低检测概率下也能够稳健地跟踪上所有目标，准确地估计目标数量和目标状态，具有很好的跟踪性能；
2.本发明提出的FMS-CPHD滤波方法的OSPA距离都明显小于传统MS-CPHD滤波器；
3.本发明提出的FMS-CPHD滤波方法在所有环境下的目标数估计都比传统MS-CPHD滤波器更加准确，尤其是当量测误差较大或检测概率较低时，说明FMS-CPHD滤波器跟踪性能有显著的提高，在较差条件下也能维持稳定的性能；
4.本发明提出的FMS-CPHD滤波方法的耗时远少于传统MS-CPHD滤波器，且运算耗时比较稳定，运算效率很高，具有很好的工程应用前景。
附图说明
图1是本发明的一种快速多传感器集势概率假设密度滤波方法的流程图；
图2是某航迹用各传感器量测组成量测分组的过程示例；
图3是用各航迹的量测分组组成全局量测划分的过程示例；
图4是本发明仿真实验的监视场景和目标真实轨迹示意图；
图5是本发明仿真实验的跟踪效果图；
图6是本发明仿真实验中FMS-CPHD与传统MS-CPHD滤波器的OSPA距离对比图；
图7是本发明仿真实验中FMS-CPHD与传统MS-CPHD滤波器的目标数估计结果对比图；
图8是本发明仿真实验中FMS-CPHD与传统MS-CPHD滤波器的算法耗时对比图。
具体实施方式
结合图1，本发明的一种快速多传感器集势概率假设密度滤波方法，包括以下步骤：
S1，场景中总共有Υ个传感器，传感器编号γ∈{1，2，...，Υ}；设初始时刻k＝0，初始化势分布(目标数分布)ρ0(n)和概率假设密度D0(x)；
S2，已有上一时刻(k-1时刻)的概率假设密度Dk-1(x)和势分布ρk-1(n)，对当前时刻(k时刻)进行预测，得到预测概率假设密度Dk|k-1(x)和预测势分布ρk|k-1(n)；
S3，获得传感器γ当前时刻的量测集合
Figure BDA0002905571730000091
其中
Figure BDA0002905571730000092
是量测状态，
Figure BDA0002905571730000093
是量测数量，通过各传感器的量测集合，计算各传感器的权重矩阵和代价矩阵；
S4，利用S3获得的代价矩阵
Figure BDA0002905571730000094
通过航迹对各传感器的量测进行分组，然后用各航迹的量测分组组成全局量测划分；
S5，通过S2获得的预测概率假设密度Dk|k-1(x)和预测势分布ρk|k-1(n)，以及S4获得的量测划分集合Ψ进行更新，获得更新的概率假设密度Dk(x)和更新势分布ρk(n)；S6，对高斯分量进行剪枝与合并，估计目标数和目标状态；
S6，对高斯分量进行剪枝与合并，估计目标数和目标状态；
S7，重复S2～S6，直到不再需要跟踪为止。
本发明的效果还可以通过以下仿真实验进一步说明：
仿真实验环境为英特尔i7 3.6Hz主频的8核CPU处理器，程序使用Matlab语言编写。
仿真实验中，提出的快速多传感器集势概率假设密度(Fast Multisensor-CPHD，FMS-CPHD)滤波器与传统多传感器集势概率假设密度(Multisensor-CPHD，MS-CPHD)滤波器进行对比。
1、仿真条件
选择一个大小为[-1000，1000]×[-1000，1000](m)的正方形二维区域作为监视区域，有三个传感器同时进行监视；
仿真过程中总共出现4个目标，目标的轨迹如图4中所示，目标在初始时刻从四个位置产生，并运动到终止时刻，目标起始位置和终止位置已在图中标识；
目标的状态表示为
Figure BDA0002905571730000101
其中[x，y]T是目标横向和纵向的位置，
Figure BDA0002905571730000102
是是目标横向和纵向的速度；
采样间隔τ＝1s，状态转移矩阵F为：
Figure BDA0002905571730000103
过程噪声协方差矩阵Q＝Gv，其中
Figure BDA0002905571730000104
且v1＝v2＝3m/s2，G由下式给出：
Figure BDA0002905571730000105
初始概率假设密度为：
Figure BDA0002905571730000106
其中初始权重
Figure BDA0002905571730000107
初始协方差
Figure BDA0002905571730000108
初始标签
Figure BDA0002905571730000109
Figure BDA00029055717300001010
是初始均值，取值分别为
Figure BDA00029055717300001011
Figure BDA00029055717300001012
和
Figure BDA00029055717300001013
之后每个时刻的新目标概率假设密度DB，k|k-1(x)与D0(x)相同；
初始势分布ρ0(n)取二项分布，最大可能目标数Nmax＝20，目标出现概率r＝0.1；
三个传感器的量测矩阵均为：
Figure BDA00029055717300001014
量测噪声协方差矩阵R＝diag([q1 2，q2 2]T)，其中量测误差取q1＝q2＝5m和q1＝q2＝20m两种情况；
杂波强度取κk(z)＝1×10-5(m2)-1，杂波平均分布在监测区域中；
三个传感器具有相同的检测概率，即
Figure BDA00029055717300001015
且检测概率取0.9和0.5两种情况；
目标存活概率为pS＝0.99。
2、仿真结果分析
图5给出了检测概率0.5，量测误差q1＝q2＝5m时，FMS-CPHD滤波器的跟踪结果。图中用不同形状的记号展示了不同传感器的量测情况。可以看出，FMS-CPHD滤波器在低检测概率下也能够稳健地跟踪上所有目标，准确地估计目标数量和目标状态，说明本发明具有很好的跟踪性能。
图6给出了FMS-CPHD滤波器和传统MS-CPHD滤波器在不同环境下进行100次蒙特卡洛仿真的平均结果，不同环境包括：检测概率0.9、测量误差q1＝q2＝5m；检测概率0.9、测量误差q1＝q2＝20m和检测概率0.5、测量误差q1＝q2＝5m三种情况。这里使用100次蒙特卡洛的平均OSPA距离作为性能评价准则。从图6中可以看出，在所有环境下，FMS-CPHD滤波器的OSPA距离都明显小于传统MS-CPHD滤波器，说明FMS-CPHD滤波器跟踪性能有显著的提高。
图7给出了FMS-CPHD滤波器和传统MS-CPHD滤波器在这100次蒙特卡洛仿真中的平均目标数估计结果。FMS-CPHD滤波器在所有环境下的目标数估计都比传统MS-CPHD滤波器更加准确，尤其是当量测误差较大或检测概率较低时，说明FMS-CPHD滤波器跟踪性能有显著的提高，在较差条件下也能维持稳定的性能。
图8给出了FMS-CPHD滤波器和传统MS-CPHD滤波器在这100次蒙特卡洛仿真中的平均算法耗时。可以看出FMS-CPHD滤波器的耗时远少于传统MS-CPHD滤波器，且运算耗时比较稳定，说明本发明运算效率很高，具有很好的工程应用前景。

Claims (12)

1.一种快速多传感器集势概率假设密度滤波方法，其特征在于，该方法包括以下步骤：
先给出统一的符号定义：(·)k表示k时刻的值，则(·)0表示初始值，(·)k|k-1表示k-1时刻对k时刻的预测值，(·)(i)表示与索引号为i的高斯分量有关的物理量，(·)γ表示与传感器γ有关的物理量，(·)l表示与目标航迹l有关的物理量；
再给出高斯分量的定义：
Figure FDA0003302110050000011
表示k时刻索引号为i的高斯分量，其权重为
Figure FDA0003302110050000012
均值为
Figure FDA0003302110050000013
协方差为
Figure FDA0003302110050000014
被标签
Figure FDA0003302110050000015
标记；索引号i∈{1,2,...,Ik}，其中Ik表示k时刻高斯分量的数量；标签
Figure FDA0003302110050000016
其中
Figure FDA0003302110050000017
表示k时刻的标签数，被同一个标签标记的高斯分量认为来自同一个目标航迹；
S1，场景中总共有γ个传感器，传感器编号γ∈{1,2,...,γ}；设初始时刻k＝0，初始化势分布ρ0(n)和概率假设密度D0(x)，具体为：
初始化势分布为ρ0(n)，其中n表示目标数为n且n∈{0,1,...,Nmax}，Nmax为最大可能目标数；ρ0(n)可以根据情况选择合适分布，这里选择使用二项分布：
Figure FDA0003302110050000018
其中C表示组合数符号，
Figure FDA0003302110050000019
r为目标出现概率；
初始化概率假设密度D0(x)的形式为：
Figure FDA00033021100500000110
S2，已有上一时刻的概率假设密度Dk-1(x)和势分布ρk-1(n)，对当前时刻进行预测，得到预测概率假设密度Dk|k-1(x)和预测势分布ρk|k-1(n)；包括以下步骤：
上一时刻的概率假设密度Dk-1(x)的形式为：
Figure FDA00033021100500000111
对当前时刻的势分布进行预测为：
Figure FDA00033021100500000112
其中ρB(·)是新目标的势分布，pS是目标存活概率，s和t表示两个整数；
对当前时刻的概率假设密度进行预测为：
Dk|k-1(x)＝DS,k|k-1(x)+DB,k|k-1(x)
其中DS,k|k-1(x)表示存活目标的概率假设密度，DB,k|k-1(x)表示新目标概率假设密度；
DS,k|k-1(x)由下式给出：
Figure FDA0003302110050000021
其中
Figure FDA0003302110050000022
是存活目标标签且
Figure FDA0003302110050000023
和
Figure FDA0003302110050000024
分别是存活目标高斯分量的均值和协方差，由下两式分别计算：
Figure FDA0003302110050000025
Figure FDA0003302110050000026
其中F和Q分别表示状态转移矩阵和过程噪声协方差矩阵，T表示矩阵转置；
DB,k|k-1(x)由下式给出：
Figure FDA0003302110050000027
其中
Figure FDA0003302110050000028
和
Figure FDA0003302110050000029
分别是新目标高斯分量的权重，均值和协方差，IB,k|k-1是新目标高斯分量数，新目标标签
Figure FDA00033021100500000210
Figure FDA00033021100500000211
是新目标数；
则预测概率假设密度Dk|k-1(x)为：
Figure FDA00033021100500000212
其中Ik|k-1＝Ik-1+IB,k|k-1，标签
Figure FDA00033021100500000213
S3，获得传感器γ当前时刻的量测集合
Figure FDA00033021100500000214
其中
Figure FDA00033021100500000215
是量测状态，
Figure FDA00033021100500000216
是量测数量，通过各传感器的量测集合，计算各传感器的权重矩阵和代价矩阵；包括以下步骤：
S3.1传感器γ的代价矩阵用
Figure FDA00033021100500000217
表示，矩阵第l行，第j列的元素
Figure FDA00033021100500000218
表示航迹l与量测
Figure FDA00033021100500000219
的关联代价，其中
Figure FDA00033021100500000220
由下式计算：
Figure FDA00033021100500000221
其中：
Figure FDA00033021100500000222
表示标记为l的高斯分量索引号的集合；
Figure FDA00033021100500000223
H是量测矩阵，R是量测噪声协方差矩阵；
Figure FDA00033021100500000224
表示传感器γ的检测概率；
κ(·)表示杂波的强度函数；
S3.2传感器γ的权重矩阵由
Figure FDA0003302110050000031
表示，其大小与
Figure FDA0003302110050000032
相同，第l行，第j列的元素
Figure FDA0003302110050000033
由下式计算：
Figure FDA0003302110050000034
S3.3计算量测
Figure FDA0003302110050000035
的更新权重：
Figure FDA0003302110050000036
设置量测剔除阈值θz＝10-3，如果
Figure FDA0003302110050000037
则认为
Figure FDA0003302110050000038
是杂波，将其所在列从代价矩阵
Figure FDA0003302110050000039
中剔除；剔除后矩阵
Figure FDA00033021100500000310
剩余
Figure FDA00033021100500000311
列；
S3.4计算传感器γ代价矩阵的行列差：
Figure FDA00033021100500000312
从各传感器代价矩阵的行列差
Figure FDA00033021100500000313
中取最大值为
Figure FDA00033021100500000314
如果Δk＞0，则将Δk+1个代价矩阵
Figure FDA00033021100500000315
纵向拼接为新的代价矩阵：
Figure FDA00033021100500000316
S3.5代价矩阵
Figure FDA00033021100500000317
的关联矩阵用
Figure FDA00033021100500000318
表示，其大小与
Figure FDA00033021100500000319
相同，其中第l行，第j列的元素为al,j，al,j等于0或1，且每行每列最多只有1个al,j等于1；则关联矩阵
Figure FDA00033021100500000320
的代价为
Figure FDA00033021100500000321
其中
Figure FDA00033021100500000322
表示此时
Figure FDA00033021100500000323
的行数，
Figure FDA00033021100500000324
表示
Figure FDA00033021100500000325
中第l行，第j列的元素；
通过寻优算法寻找代价最小的关联矩阵
Figure FDA00033021100500000326
然后求矩阵
Figure FDA00033021100500000327
与
Figure FDA00033021100500000328
的Hadamard积，得到最终的代价矩阵
Figure FDA00033021100500000329
其中符号*表示Hadamard积运算；
S4，利用S3获得的代价矩阵
Figure FDA00033021100500000330
通过航迹对各传感器的量测进行分组，然后用各航迹的量测分组组成全局量测划分，包括以下步骤：
S4.1用每个航迹对各传感器的量测进行分组，具体过程为：
对于航迹l，依次从传感器γ∈{1,2,...,γ}的量测集合
Figure FDA00033021100500000331
中选择一个量测，且只有在传感器γ的代价矩阵
Figure FDA00033021100500000332
中
Figure FDA00033021100500000333
的量测可以被选择，不选择量测则表示传感器γ漏检该目标航迹l；完成后将所有选择的量测组成航迹l的第1个量测分组
Figure FDA00033021100500000334
S4.2重复S4.1的过程，找到航迹l所有可能的量测分组
Figure FDA0003302110050000041
其中Ll是航迹l的量测分组数；
S4.3对所有的航迹进行S4.1～S4.2的过程，找到所有航迹的量测分组；
S4.4然后用各航迹的量测分组组成全局量测划分，具体过程为：
依次从航迹
Figure FDA0003302110050000042
的量测分组中选择一个分组
Figure FDA0003302110050000043
其中ll∈{1,2,...,Ll}，将这些量测分组组成一个全局量测划分
Figure FDA0003302110050000044
其中
Figure FDA0003302110050000045
表示来自所有传感器的杂波组成的集合；用
Figure FDA0003302110050000046
表示量测分组
Figure FDA0003302110050000047
中来自传感器γ的量测数量，则ψ中来自传感器γ的量测数量表示为
Figure FDA0003302110050000048
ψ中总的量测数量表示为
Figure FDA0003302110050000049
S4.5重复该过程，找到所有可能的量测划分，所有量测划分组成集合Ψ；
S5，通过S2获得的预测概率假设密度Dk|k-1(x)和预测势分布ρk|k-1(n)，以及S4获得的量测划分集合Ψ进行更新，获得更新的概率假设密度Dk(x)和更新势分布ρk(n)；包括以下步骤：
S5.1更新的势分布ρk(n)由下式计算：
Figure FDA00033021100500000410
其中：
Figure FDA00033021100500000411
Figure FDA00033021100500000412
是传感器γ的杂波的势分布；
Figure FDA00033021100500000413
Figure FDA00033021100500000414
Figure FDA00033021100500000415
P表示排列数符号，
Figure FDA00033021100500000416
Figure FDA00033021100500000417
表示量测分组
Figure FDA00033021100500000418
的似然比，由下式计算：
Figure FDA0003302110050000051
其中：
Figure FDA0003302110050000052
cγ(·)是杂波的空间分布函数；
S5.2索引号为i的高斯分量在量测分组
Figure FDA0003302110050000053
中的归一化伪似然
Figure FDA0003302110050000054
由下式计算：
Figure FDA0003302110050000055
S5.3更新的概率假设密度Dk(x)为：
Dk(x)＝DE,k(x)+DU,k(x)
其中DE,k(x)和DU,k(x)分别表示遗留概率假设密度和量测更新概率假设密度，分别为：
Figure FDA0003302110050000056
Figure FDA0003302110050000057
其中：
Figure FDA0003302110050000058
而α0由下式计算：
Figure FDA0003302110050000059
Figure FDA00033021100500000510
表示全局量测划分ψ中量测分组
Figure FDA00033021100500000511
更新的概率假设密度，由下式给出：
Figure FDA00033021100500000512
其中比例系数
Figure FDA00033021100500000513
由下式计算：
Figure FDA0003302110050000061
量测分组
Figure FDA0003302110050000062
更新的高斯分量权重
Figure FDA0003302110050000063
由下式计算：
Figure FDA0003302110050000064
其中αψ由下式计算：
Figure FDA0003302110050000065
Figure FDA0003302110050000066
是量测
Figure FDA0003302110050000067
更新的高斯分量均值，
Figure FDA0003302110050000068
是更新协方差，分别由下式计算：
Figure FDA0003302110050000069
Figure FDA00033021100500000610
Figure FDA00033021100500000611
则量测更新概率假设密度DU,k(x)为：
Figure FDA00033021100500000612
其中
Figure FDA00033021100500000613
新标签
Figure FDA00033021100500000614
由量测分组
Figure FDA00033021100500000615
决定，
Figure FDA00033021100500000616
是航迹l选择的量测分组，则
Figure FDA00033021100500000617
(5.4)重新整理DU,k(x)中高斯分量的索引号，可将DU,k(x)整理为：
Figure FDA00033021100500000618
其中
Figure FDA00033021100500000619
为所有ψ中量测数量之和；
将DE,k(x)重新整理为：
Figure FDA00033021100500000620
其中
Figure FDA00033021100500000621
最终，所有高斯分量组成当前时刻更新的概率假设密度Dk(x)：
Figure FDA0003302110050000071
其中Ik＝Ik|k-1+Ik|k-1·|Ψ|；
S6，对高斯分量进行剪枝与合并，估计目标数和目标状态，包括以下步骤：
S6.1设置高斯分量剔除阈值wP＝10-5，如果
Figure FDA0003302110050000072
将权重
Figure FDA0003302110050000073
对应的高斯分量删除；
S6.2如果具有相同标签的高斯分量之间的距离小于距离阈值dM＝4m，将这些高斯分量合并；
S6.3目标数估计
Figure FDA0003302110050000074
就是使更新势分布取最大值的整数，即：
Figure FDA0003302110050000075
S6.4航迹l的高斯分量权重和为：
Figure FDA0003302110050000076
将航迹高斯分量权重和按从大到小排序，取前
Figure FDA0003302110050000077
条航迹，然后从每个航迹对应的高斯分量中选择权重最大的高斯分量，其均值就是目标状态的估计；
S7，重复S2～S6，直到不再需要跟踪为止。
2.一种根据权利要求1所述快速多传感器集势概率假设密度滤波方法，其特征在于：S3.3中，量测剔除阈值θz设置为不大于10-3的数。
3.一种根据权利要求2所述快速多传感器集势概率假设密度滤波方法，其特征在于：S3.3中，量测剔除阈值θz取10-3。
4.一种根据权利要求1所述快速多传感器集势概率假设密度滤波方法，其特征在于：S3.5中，寻优算法为匈牙利寻优算法。
5.一种根据权利要求1所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.1中，高斯分量剔除阈值wP设置为不大于10-5的数。
6.一种根据权利要求5所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.1中，高斯分量剔除阈值wP取10-5。
7.一种根据权利要求1所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.2中，合并距离阈值dM根据跟踪场景的大小决定。
8.一种根据权利要求7所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.2中，合并距离阈值dM取4m。
9.一种根据权利要求1所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.3中，高斯分量总数阈值Imax设置为不小于100的整数。
10.一种根据权利要求9所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.3中，高斯分量总数阈值Imax取200。
11.一种根据权利要求1所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.5中，航迹剔除阈值θl设置为不大于10-3的数。
12.一种根据权利要求11所述快速多传感器集势概率假设密度滤波方法，其特征在于：S6.5中，航迹剔除阈值θl取10-4。
CN202110069747.2A 2021-01-19 2021-01-19 一种快速多传感器集势概率假设密度滤波方法 Active CN112906743B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110069747.2A CN112906743B (zh) 2021-01-19 2021-01-19 一种快速多传感器集势概率假设密度滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110069747.2A CN112906743B (zh) 2021-01-19 2021-01-19 一种快速多传感器集势概率假设密度滤波方法

Publications (2)

Publication Number Publication Date
CN112906743A CN112906743A (zh) 2021-06-04
CN112906743B true CN112906743B (zh) 2021-11-19

Family

ID=76115774

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110069747.2A Active CN112906743B (zh) 2021-01-19 2021-01-19 一种快速多传感器集势概率假设密度滤波方法

Country Status (1)

Country Link
CN (1) CN112906743B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103679753A (zh) * 2013-12-16 2014-03-26 深圳大学 一种概率假设密度滤波器的轨迹标识方法及轨迹标识系统
CN106021194A (zh) * 2016-05-19 2016-10-12 哈尔滨工业大学 一种多传感器多目标跟踪偏差估计方法
CN106405538A (zh) * 2016-09-13 2017-02-15 深圳大学 一种适用于杂波环境的多目标跟踪方法及跟踪系统
CN110376547A (zh) * 2019-07-20 2019-10-25 中国人民解放军国防科技大学 基于二阶统计量的近场源定位方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2010315600A1 (en) * 2009-10-26 2012-05-10 Abbott Laboratories Diagnostic methods for determining prognosis of non-small cell lung cancer

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103679753A (zh) * 2013-12-16 2014-03-26 深圳大学 一种概率假设密度滤波器的轨迹标识方法及轨迹标识系统
CN106021194A (zh) * 2016-05-19 2016-10-12 哈尔滨工业大学 一种多传感器多目标跟踪偏差估计方法
CN106405538A (zh) * 2016-09-13 2017-02-15 深圳大学 一种适用于杂波环境的多目标跟踪方法及跟踪系统
CN110376547A (zh) * 2019-07-20 2019-10-25 中国人民解放军国防科技大学 基于二阶统计量的近场源定位方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
An_improved_Multitarget_Multi-Bernoulli_filter_with_cardinality_corrected;Zhejun Lu;Weidong Hu;Hua Yu;Thia Kirubarajan;《Institute of Electrical and Electronics Engineers》;IEEE;20160708;全文 *
基于概率假设密度滤波的多扩展目标跟踪技术;彭聪;王杰贵;张坤;《电子信息对抗技术》;20180630;全文 *
空间碎片群目标状态估计理论与方法研究;卢哲俊;《CNKI优秀硕士学位论文全文库》;CNKI优秀硕士学位论文全文库;20171031;全文 *

Also Published As

Publication number Publication date
CN112906743A (zh) 2021-06-04

Similar Documents

Publication Publication Date Title
CN110503071B (zh) 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法
CN107831490A (zh) 一种改进的多扩展目标跟踪方法
CN111722214B (zh) 雷达多目标跟踪phd实现方法
CN106991691A (zh) 一种适用于摄像机网络下的分布式目标跟踪方法
CN107063259A (zh) 一种航迹关联方法及电子设备
CN106934324A (zh) 基于简化多假设算法的雷达数据关联方法
CN103985120A (zh) 一种遥感图像多目标关联的方法
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
CN109214432B (zh) 一种多传感器多目标联合检测、跟踪与分类方法
Chavali et al. Hierarchical particle filtering for multi-modal data fusion with application to multiple-target tracking
CN109509207B (zh) 一种对点目标和扩展目标进行无缝跟踪的方法
CN113344970B (zh) 基于多伯努利的非规则多扩展目标联合跟踪与分类方法
CN111274529B (zh) 一种鲁棒的高斯逆威沙特phd多扩展目标跟踪算法
CN114002667A (zh) 一种基于随机矩阵法的多近邻扩展目标跟踪算法
CN109671096B (zh) 一种时空近邻目标检测及网格聚类量测划分下的多扩展目标跟踪方法
CN112906743B (zh) 一种快速多传感器集势概率假设密度滤波方法
CN111504327B (zh) 一种基于航迹平滑技术的广义标签多伯努利多目标跟踪方法
CN110045363B (zh) 基于相对熵的多雷达航迹关联方法
CN116047495B (zh) 一种用于三坐标雷达的状态变换融合滤波跟踪方法
CN116224320B (zh) 一种极坐标系下处理多普勒量测的雷达目标跟踪方法
CN116736290A (zh) 一种基于空中生物运动态势特征辅助的多目标跟踪方法
CN104467742A (zh) 基于高斯混合模型的传感器网络分布式一致性粒子滤波器
Yang et al. Multiple extended target tracking algorithm based on Gaussian surface matrix
CN111104985B (zh) 一种异步航迹关联的加权滑窗方法
CN111811515A (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