CN108732564B - 一种双雷达修正序贯高斯混合概率假设密度滤波方法 - Google Patents

一种双雷达修正序贯高斯混合概率假设密度滤波方法 Download PDF

Info

Publication number
CN108732564B
CN108732564B CN201810203353.XA CN201810203353A CN108732564B CN 108732564 B CN108732564 B CN 108732564B CN 201810203353 A CN201810203353 A CN 201810203353A CN 108732564 B CN108732564 B CN 108732564B
Authority
CN
China
Prior art keywords
gaussian
radar
gaussian component
measurement
ith
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
CN201810203353.XA
Other languages
English (en)
Other versions
CN108732564A (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.)
CETC 28 Research Institute
Original Assignee
CETC 28 Research Institute
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 CETC 28 Research Institute filed Critical CETC 28 Research Institute
Priority to CN201810203353.XA priority Critical patent/CN108732564B/zh
Publication of CN108732564A publication Critical patent/CN108732564A/zh
Application granted granted Critical
Publication of CN108732564B publication Critical patent/CN108732564B/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
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00

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

本发明公开了一种双雷达修正序贯高斯混合概率假设密度滤波方法，传统双雷达序贯GM‑PHD仅适用于测量目标位于两个雷达的公共测量区域的情况，当目标未处于上述区域时，可能会出现目标丢失问题。本发明基于有限统计学理论，通过对各个雷达测量值对应的高斯成分进行预测、更新、修剪融合、维持融合、目标状态提取，实现多目标跟踪，且不会丢失处于非公共测量区域的目标，拓宽了序贯GM‑PHD的应用范围。与传统方法相比，计算复杂性变化不大。

Description

一种双雷达修正序贯高斯混合概率假设密度滤波方法
技术领域
本发明涉及多传感器多目标跟踪技术领域，尤其涉及一种双雷达修正序贯高斯混合概率假设密度滤波方法。
背景技术
20世纪50年代，随着雷达应用环境的复杂化，要求雷达能够同时跟踪多个目标，多目标跟踪的概念随之提出。经过几十年的研究，多目标跟踪技术理论发展迅速，越来越多的优良算法被提出，并广泛应用于军事和民事的各个领域，如军事情报收集、敌情预警、工业过程监控、空中交通管制等。
工程中常用的多目标跟踪算法主要有最近邻法(Nearest Neighbor,NN)、联合概率数据互联法(Joint Probabilistic Data Association，JPDA)以及多假设跟踪算法(Multiple Hypothesis Tracking,MHT)等。这些算法都需要通过数据关联技术将多目标跟踪问题转化为单目标跟踪问题，实现量测数据与已有航迹的准确关联，但数据关联本身是NP难的，也就是说引入数据关联大大增加了多目标跟踪问题的复杂度。因此，需要寻找一种不需要进行数据关联的多目标跟踪算法。
基于有限集统计(Finite Set Statistics,FISST)理论，一种不需要进行数据关联的多目标跟踪算法--概率假设密度(Probability Hypothesis Density,PHD)滤波--被提出，为解决多目标跟踪提供了一种新方法。由于不需要进行数据关联，PHD滤波的计算复杂度比传统的多目标跟踪算法要小，但其迭代式中包含多个一般情况下没有封闭形式解的积分，不便于工程应用。为此，基于线性高斯假设，澳大利亚科学家提出了一种高斯混合概率假设密度(Gaussian Mixture Probability Hypothesis Density,GM-PHD)滤波，使PHD滤波的工程应用向前迈出了一大步。
随着对跟踪目标实时性和准确性的要求越来越高，多传感器信息融合迅速发展。上面提到的各种滤波算法都是针对单传感器多目标跟踪，不能直接用于多传感器多目标跟踪。为解决多传感器多目标跟踪问题，许多人对传统多目标跟踪算法进行拓展，使其适用于多传感器目标跟踪。序贯GM-PHD就是一种基于GM-PHD的多传感器多目标跟踪滤波算法。当目标仅被部分传感器测量到时，序贯GM-PHD会出现目标丢失问题，这就大大限制了其使用范围。
实际系统中，两个雷达的测量范围一般不会完全重合，而目标可能随机从任意地点出现，当目标从一个雷达测量范围的边缘出现时，传统序贯GM-PHD将会丢失目标，严重影响目标跟踪的准确性。在战场上，不能及时发现目标，可能会导致误判，影响战斗的结果。
发明内容
为了解决双雷达系统中目标位于非公共测量区域时，序贯GM-PHD滤波会导致目标丢失的问题，本发明提供一种保持仅被单雷达测量到的目标的方法。
现有的序贯GM-PHD滤波只能针对目标同时被两个雷达测量到的情况，仅被一个雷达测量到的目标会丢失。为了使序贯GM-PHD滤波的应用可以不受目标位置的限制，本发明在传统序贯GM-PHD滤波中加入维持融合操作，从而实现了双雷达系统中的多目标正确跟踪。
本发明所采用的技术方案是：双雷达系统中修正序贯GM-PHD滤波。本发明假设雷达测量时间同步，上一时刻的目标状态强度的高斯成分已知，序贯算法先使用雷达1的测量数据，后使用雷达2的测量数据，具体包括如下步骤：
步骤1，初始化：雷达1和雷达2组成雷达系统，假设雷达1和雷达2测量时间同步，建立系统方程和测量方程，通过传感器获得目标初始时刻的状态，包含目标的位置和速度信息，得到第0个周期的目标状态强度的高斯成分集合；
步骤2，雷达1高斯成分预测：以k-1表示雷达系统前一扫描周期，以k表示雷达当前扫描周期，记雷达扫描周期为T，利用系统方程将k-1周期目标状态强度的高斯成分集合进行预测，得到k周期雷达1的预测高斯成分集合，该预测高斯成分集合中包括存活高斯成分预测项、新出生高斯成分项和分裂高斯成分项；
步骤3，雷达1高斯成分测量更新；
步骤4，雷达1高斯成分融合修剪：删除测量更新后的高斯成分集合中权重充分小的高斯项，将距离充分小的高斯成分合并成一个；
步骤5，雷达2高斯成分测量更新；
步骤6，雷达2高斯成分融合修剪：删除测量更新后的高斯成分集合中权重充分小的高斯项，将距离充分小的高斯成分合并成一个；
步骤7，高斯成分维持融合：对于步骤4得到的雷达1的修剪融合后的高斯成分集合中的每一个元素，从步骤6得到的雷达2修剪融合后的高斯成分集合中找出与其距离最近的元素，进行维持融合，得到维持融合后的高斯成分集合；
步骤8，目标状态提取，将步骤7得到的高斯成分集合中权重大于0.5的高斯成份均值提取出来作为当前时刻的目标状态估计值。
步骤1包括如下步骤：
步骤1-1，建立如下系统方程：
Xk＝Fk*Xk-1+wk-1 (公式一)
其中，
Figure BDA0001595166130000031
表示k时刻目标状态向量，xk为x方向的位置，yk为y方向的位置，
Figure BDA0001595166130000032
为x方向的速度，
Figure BDA0001595166130000033
为y方向的速度，Fk为k时刻状态转移矩阵，wk为k时刻系统噪声，wk服从均值为0，协方差为Qk的正态分布，Qk为k时刻正态分布的协方差；
步骤1-2，为生成雷达1和雷达2的测量数据，建立如下两个测量方程：
Figure BDA0001595166130000034
其中，
Figure BDA0001595166130000035
为雷达1的测量值，
Figure BDA0001595166130000036
为雷达2的测量值，两个雷达的测量矩阵满足
Figure BDA0001595166130000037
雷达1和雷达2的测量噪声分别为具有协方差
Figure BDA0001595166130000038
的零均值、白色高斯测量噪声
Figure BDA0001595166130000039
和具有协方差
Figure BDA00015951661300000310
的零均值、白色高斯测量噪声
Figure BDA00015951661300000311
分别表示为
Figure BDA00015951661300000312
N表示正态分布，
Figure BDA00015951661300000313
Figure BDA00015951661300000318
分别为雷达1的测量噪声的标准差和雷达2的测量噪声的标准差，I2为2阶单位矩阵；
步骤1-3，通过雷达1获得目标初始时刻的状态，包含目标的位置和速度信息，得到第0个周期的目标状态强度的高斯成分集合为
Figure BDA00015951661300000314
其中，
Figure BDA00015951661300000315
为0时刻第i个高斯成分的权重，
Figure BDA00015951661300000316
为0时刻第i个高斯成分的状态向量，
Figure BDA00015951661300000317
为0时刻第i个高斯成分的方差阵，J0为0时刻高斯成分的数量。
步骤1-1中，状态转移矩阵Fk和噪声方差阵Qk分别如下：
Figure BDA0001595166130000041
其中，T为测量周期，
Figure BDA0001595166130000042
为系统噪声的标准差，F、Q为中间变量。
步骤2中，得到k时刻雷达1的预测高斯成分集合为
Figure BDA0001595166130000043
其中，
Figure BDA0001595166130000044
为k时刻第i个预测高斯成分的权重，
Figure BDA0001595166130000045
为k时刻第i个预测高斯成分的状态向量，
Figure BDA0001595166130000046
为k时刻第i个预测高斯成分的方差阵，Jk|k-1为k时刻预测高斯成分的数量。
步骤3包括如下步骤：
步骤3-1，根据雷达1的预测高斯成分集合计算预测测量值
Figure BDA0001595166130000047
预测测量误差协方差阵
Figure BDA0001595166130000048
增益矩阵
Figure BDA0001595166130000049
和估计误差协方差阵
Figure BDA00015951661300000410
公式如下：
Figure BDA00015951661300000411
其中，I为单位矩阵；
步骤3-2，对k时刻雷达1的每一个测量值
Figure BDA00015951661300000412
对预测高斯成分进行测量更新得到多目标状态集合矩
Figure BDA00015951661300000413
表示为：
Figure BDA00015951661300000414
其中，i为下标，
Figure BDA00015951661300000415
表示将第i个预测高斯成分作为以x为变量，以
Figure BDA00015951661300000416
为均值、以
Figure BDA00015951661300000417
为方差阵的高斯分布，
Figure BDA00015951661300000418
为高斯分布的权重，
Figure BDA00015951661300000419
表示用第j个测量值更新第i个预测高斯成分后得到的更新后的以x为变量、以
Figure BDA00015951661300000420
为均值、以
Figure BDA00015951661300000421
为方差阵的高斯分布，其权重为
Figure BDA00015951661300000422
Figure BDA0001595166130000051
其中，
Figure BDA0001595166130000052
为第i个高斯成分的状态，
Figure BDA0001595166130000053
为k时刻第l个预测高斯成分的权重，
Figure BDA0001595166130000054
为第i个高斯成分的方差阵，kk(zj)为第j个测量值的杂波强度函数，
Figure BDA0001595166130000055
为雷达1由第1至k测量周期所有测量构成的集合，
Figure BDA0001595166130000056
为雷达1检测到目标的概率；
步骤3-3，得到测量更新后的高斯成分集合
Figure BDA0001595166130000057
其中，
Figure BDA0001595166130000058
为k时刻第i个更新后高斯成分的权重，
Figure BDA0001595166130000059
为k时刻第i个更新后高斯成分的状态向量，
Figure BDA00015951661300000510
为k时刻第i个更新后高斯成分的方差阵，Jk|k为k时刻更新后高斯成分的数量。
步骤4包括如下步骤：
步骤4-1，删除测量更新后的高斯成分集合
Figure BDA00015951661300000511
中权重充分小的高斯项，即删除权重
Figure BDA00015951661300000512
的高斯项，其中
Figure BDA00015951661300000522
为裁减门限，一般可取为0.00001；
步骤4-2，将距离充分小的高斯成分合并成一个，即将
Figure BDA00015951661300000513
的高斯项合并成一个，其中
Figure BDA00015951661300000514
为合并门限，一般可取为4，距离di,j定义为
Figure BDA00015951661300000515
高斯成分合并方法如下：
Figure BDA00015951661300000516
其中L为所有可合并高斯成分的指标集合，
Figure BDA00015951661300000517
为k时刻第i个合并后高斯成分的权重，
Figure BDA00015951661300000518
为k时刻第i个合并后高斯成分的状态向量，
Figure BDA00015951661300000519
为k时刻第i个合并后高斯成分的方差阵，j为k时刻L中可合并高斯成分的下标；
步骤4-3，得到雷达1修剪融合后的高斯成分集合为
Figure BDA00015951661300000520
其中，
Figure BDA00015951661300000521
为修剪融合后的第i个高斯成分，
Figure BDA0001595166130000061
为其状态，
Figure BDA0001595166130000062
为其协方差阵，Jk,1为雷达1修剪融合后高斯成分的个数。为进行序贯操作，将该高斯成分集合作为雷达2的预测高斯成分集合。
步骤5包括如下步骤：
步骤5-1，根据雷达2的预测高斯成分集合计算预测测量值
Figure BDA0001595166130000063
预测测量误差协方差阵
Figure BDA0001595166130000064
增益矩阵
Figure BDA0001595166130000065
估计误差协方差阵
Figure BDA0001595166130000066
公式如下：
Figure BDA0001595166130000067
步骤5-2，对k时刻雷达2的每一个测量值
Figure BDA0001595166130000068
对预测高斯成分进行测量更新得到多目标状态集合矩
Figure BDA0001595166130000069
表示为：
Figure BDA00015951661300000610
其中，
Figure BDA00015951661300000611
其中，kk(z)为杂波强度函数，
Figure BDA00015951661300000612
为雷达2由第1至k测量周期所有测量构成的集合，
Figure BDA00015951661300000613
为雷达2检测到目标的概率；
步骤5-3，得到测量更新后的高斯成分集合
Figure BDA00015951661300000614
步骤6包括如下步骤：
步骤6-1，删除测量更新后的高斯成分集合
Figure BDA00015951661300000615
中权重充分小的高斯项，即删除权重
Figure BDA00015951661300000616
的高斯项，其中
Figure BDA00015951661300000617
为裁减门限，一般可取为0.00001；
步骤6-2，将距离充分小的高斯成分合并成一个，即将
Figure BDA0001595166130000071
的高斯成分合并成一个，其中
Figure BDA0001595166130000072
为合并门限，一般可取为4，距离di,j定义为
Figure BDA0001595166130000073
高斯成分合并方法如下：
Figure BDA0001595166130000074
其中L为所有可合并高斯成分的指标集合；
步骤6-3，得到雷达2修剪融合后的高斯成分集合为
Figure BDA0001595166130000075
其中，
Figure BDA0001595166130000076
为修剪融合后的第i个高斯成分，
Figure BDA0001595166130000077
为其状态，
Figure BDA0001595166130000078
为其协方差阵，Jk,2为雷达2修剪融合后高斯成分的个数。
步骤7包括如下步骤：
步骤7-1，对于步骤4得到的雷达1的修剪融合后的高斯成分集合I1，即
Figure BDA0001595166130000079
中的每一个元素，从步骤6得到的雷达2修剪融合后的高斯成分集合I2，即
Figure BDA00015951661300000710
中找出与其距离最近的元素，设
Figure BDA00015951661300000711
为最近距离，M为维持融合门限，一般可取为4，若dj>M，则将I1中元素
Figure BDA00015951661300000712
添加到I2中，否则，通过两个元素中的权值项进行比较，将两个元素融合为一个元素，融合的具体过程为：
若I1中元素权值不小于0.5而I2中元素权值小于0.5，则用I1中元素替换I2中元素，若I1中元素权值小于0.5而I2中元素权值不小于0.5，则不进行操作，若两元素权值均小于或均不小于0.5，则进行加权平均，公式如下：
Figure BDA00015951661300000713
其中
Figure BDA0001595166130000081
为I1中元素，
Figure BDA0001595166130000082
为I2中元素，w1，w2和m为中间变量；
步骤7-2，得到维持融合后的高斯成分集合为
Figure BDA0001595166130000083
其中，
Figure BDA0001595166130000084
为维持融合后第i个高斯成分的权重，
Figure BDA0001595166130000085
为其状态，
Figure BDA0001595166130000086
为其协方差阵，Jk为维持融合后高斯成分的数量。
有益效果：本发明采用GM-PHD滤波，不需要进行数据关联，相比其他传统滤波算法计算量小，并且，不论目标仅被一个雷达测量到还是同时被两个雷达测量到都能够正确跟踪目标，不会出现目标丢失问题。
附图说明
下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述或其他方面的优点将会变得更加清楚。
图1是本发明的多目标跟踪方法原理框图。
图2是本发明实例目标编队飞行的轨迹图。
图3是传统方法的目标跟踪轨迹。
图4是本发明方法的目标跟踪轨迹。
图5a～图5c是本发明方法的目标跟踪轨迹局部放大图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
本发明的双雷达系统中修正序贯GM-PHD滤波能给出作为滤波器输出的高斯项，在滤波器预测、更新、裁剪合并和维持融合后，滤波器输出目标状态。本实施例中系统方程为：
Xk＝Fk*Xk-1+wk-1 (公式一)
其中
Figure BDA0001595166130000087
表示k时刻目标状态向量，每个分量对应表示目标的位置和速度，系统噪声wk～N(0,Qk)，状态转移矩阵及噪声方差阵分别为
Figure BDA0001595166130000091
测量方程分别为
Figure BDA0001595166130000092
其中
Figure BDA0001595166130000093
测量噪声
Figure BDA0001595166130000094
Figure BDA0001595166130000095
如图1所示，本发明主要包括：初始化模块、雷达1高斯成分预测模块、雷达1高斯成分测量更新模块、雷达1高斯成分融合修剪模块、雷达2高斯成分测量更新模块、雷达2高斯成分融合修剪模块、高斯成分维持融合模块、目标状态提取模块。结合流程图说明具体实现步骤为：
步骤1：初始化，通过传感器获得目标初始时刻的状态，包含目标的位置和速度信息,得到第0个周期的目标状态强度的高斯成分集合为
Figure BDA0001595166130000096
步骤2：雷达1高斯成分预测，以k-1表示雷达系统前一扫描周期，以k表示雷达当前扫描周期，记雷达扫描周期为T，利用系统方程将k-1周期目标状态强度的高斯成分集合进行预测，得到k周期预测高斯成分集合
Figure BDA0001595166130000097
该预测高斯成分集合中包括存活高斯成分预测项、新出生高斯成分项和分裂高斯成分项；
步骤3：雷达1高斯成分测量更新，首先根据雷达1的预测高斯成分集合计算预测测量值、预测测量误差协方差阵、增益矩阵、估计误差协方差阵如下:
Figure BDA0001595166130000098
对k时刻雷达1的每一个测量值
Figure BDA0001595166130000099
对预测高斯成分进行测量更新得到多目标状态集合矩，表示为：
Figure BDA0001595166130000101
其中
Figure BDA0001595166130000102
其中，kk(z)为杂波强度函数，
Figure BDA0001595166130000103
为雷达1由第1至k测量周期所有测量构成的集合，
Figure BDA0001595166130000104
为雷达1检测到目标的概率。于是得测量更新后的高斯成分集合
Figure BDA0001595166130000105
步骤4：雷达1高斯成分融合修剪，删除测量更新后的高斯成分集合
Figure BDA0001595166130000106
中权重充分小的高斯项，即删除权重
Figure BDA0001595166130000107
的高斯项，其中
Figure BDA0001595166130000108
为裁减门限，然后将距离充分小的高斯成分合并成一个，即将
Figure BDA0001595166130000109
的高斯项合并成一个，其中
Figure BDA00015951661300001010
为合并门限，距离di,j定义为
Figure BDA00015951661300001011
高斯成分合并方法如下
Figure BDA00015951661300001012
其中L为所有可合并高斯成分的指标集合。于是得到雷达1修剪融合后的高斯成分集合为
Figure BDA00015951661300001013
为进行序贯操作，将其作为雷达2的预测高斯成分集合。
步骤5：雷达2高斯成分测量更新，首先根据雷达2的预测高斯成分集合计算预测测量值、预测测量误差协方差阵、增益矩阵、估计误差协方差阵如下:
Figure BDA00015951661300001014
对k时刻雷达2的每一个测量值
Figure BDA0001595166130000111
对预测高斯成分进行测量更新得到多目标状态集合矩，表示为：
Figure BDA0001595166130000112
其中
Figure BDA0001595166130000113
其中，kk(z)为杂波强度函数，
Figure BDA0001595166130000114
为雷达2由第1至k测量周期所有测量构成的集合，
Figure BDA0001595166130000115
为雷达2检测到目标的概率。于是得测量更新后的高斯成分集合
Figure BDA0001595166130000116
步骤6：雷达2高斯成分融合修剪，删除测量更新后的高斯成分集合
Figure BDA0001595166130000117
中权重充分小的高斯项，即删除权重
Figure BDA0001595166130000118
的高斯项，其中
Figure BDA0001595166130000119
为裁减门限，然后将距离充分小的高斯成分合并成一个，即将
Figure BDA00015951661300001110
的高斯成分合并成一个，其中
Figure BDA00015951661300001111
为合并门限，距离di,j定义为
Figure BDA00015951661300001112
高斯成分合并方法如下
Figure BDA00015951661300001113
其中L为所有可合并高斯成分的指标集合。于是得到雷达2修剪融合后的高斯成分集合为
Figure BDA00015951661300001114
步骤7：高斯成分维持融合，对于步骤4得到的雷达1的修剪融合后的高斯成分集合I1，即
Figure BDA00015951661300001115
中的每一个元素，从步骤7得到的雷达2修剪融合后的高斯成分集合I2，即
Figure BDA00015951661300001116
中找出与其距离最近的元素。设
Figure BDA0001595166130000121
为最近距离，M为维持融合门限，若dj>M，则将I1中元素
Figure BDA0001595166130000122
添加到I2中，否则，通过两个元素中的权值项进行比较，将两个元素融合为一个元素。融合的具体过程为，若I1中元素权值不小于0.5而I2中元素权值小于0.5，则用I1中元素替换I2中元素，若I1中元素权值小于0.5而I2中元素权值不小于0.5，则不进行操作，若两元素权值均小于或均不小于0.5，则进行加权平均，具体为
Figure BDA0001595166130000123
其中
Figure BDA0001595166130000124
为I1中元素，
Figure BDA0001595166130000125
为I2中元素。于是得到维持融合后的高斯成分集合为
Figure BDA0001595166130000126
步骤8：目标状态提取，将权重大于0.5的高斯成份均值提取出来作为当前时刻的目标状态估计值。
在既有新目标出现和已有目标消失的情况下，对图2所示的仿真数据进行处理，传统方法得到的跟踪结果如图3所示，本发明得到的跟踪结果如图4所示，图5a为图4中并行目标运动轨迹的局部放大效果，图5b为图4中交叉目标运动轨迹的局部放大效果，图5c为图4中分裂目标运动轨迹的局部放大效果。从图4以及图5a～图5c中可以看出，本发明提出的跟踪方法能探测出观测空间的全部6批目标并能进行有效跟踪。
和现有序贯GM-PHD滤波相比，本发明的特点是：通过对两个雷达修剪融合后的高斯成分进行维持融合，可以正常跟踪从雷达探测范围内任意已知位置起始的目标，从而解决了目标位于雷达测量非公共区域时传统序贯GM-PHD会丢失目标的问题，提高了序贯GM-PHD滤波的工程使用价值。
本发明提供了一种双雷达修正序贯高斯混合概率假设密度滤波方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (2)

1.一种双雷达修正序贯高斯混合概率假设密度滤波方法，其特征在于，包括如下步骤：
步骤1，初始化：雷达1和雷达2组成雷达系统，假设雷达1和雷达2测量时间同步，建立系统方程和测量方程，通过传感器获得目标初始时刻的状态，包含目标的位置和速度信息，得到第0个周期的目标状态强度的高斯成分集合；
步骤2，雷达1高斯成分预测：以k-1表示雷达系统前一扫描周期，以k表示雷达当前扫描周期，记雷达扫描周期为T，利用系统方程将k-1周期目标状态强度的高斯成分集合进行预测，得到k周期雷达1的预测高斯成分集合，该预测高斯成分集合中包括存活高斯成分预测项、新出生高斯成分项和分裂高斯成分项；
步骤3，雷达1高斯成分测量更新；
步骤4，雷达1高斯成分融合修剪：删除测量更新后的高斯成分集合中权重充分小的高斯项，将距离充分小的高斯成分合并成一个；
步骤5，雷达2高斯成分测量更新；
步骤6，雷达2高斯成分融合修剪：删除测量更新后的高斯成分集合中权重充分小的高斯项，将距离充分小的高斯成分合并成一个；
步骤7，高斯成分维持融合：对于步骤4得到的雷达1的修剪融合后的高斯成分集合中的每一个元素，从步骤6得到的雷达2修剪融合后的高斯成分集合中找出与其距离最近的元素，进行维持融合，得到维持融合后的高斯成分集合；
步骤8，目标状态提取，将步骤7得到的高斯成分集合中权重大于0.5的高斯成份均值提取出来作为当前时刻的目标状态估计值；
步骤1包括如下步骤：
步骤1-1，建立如下系统方程：
Xk＝Fk*Xk-1+wk-1 (公式一)
其中，
Figure FDA0002341981190000011
表示k时刻目标状态向量，xk为x方向的位置，yk为y方向的位置，
Figure FDA0002341981190000012
为x方向的速度，
Figure FDA0002341981190000013
为y方向的速度，Fk为k时刻状态转移矩阵，wk为k时刻系统噪声，wk服从均值为0，协方差为Qk的正态分布，Qk为k时刻正态分布的协方差；
步骤1-2，为生成雷达1和雷达2的测量数据，建立如下两个测量方程：
Figure FDA0002341981190000021
其中，
Figure FDA0002341981190000022
为雷达1的测量值，
Figure FDA0002341981190000023
为雷达2的测量值，两个雷达的测量矩阵满足
Figure FDA0002341981190000024
雷达1和雷达2的测量噪声分别为具有协方差
Figure FDA0002341981190000025
的零均值、白色高斯测量噪声
Figure FDA0002341981190000026
和具有协方差
Figure FDA0002341981190000027
的零均值、白色高斯测量噪声
Figure FDA0002341981190000028
分别表示为
Figure FDA0002341981190000029
N表示正态分布，
Figure FDA00023419811900000210
Figure FDA00023419811900000211
分别为雷达1的测量噪声的标准差和雷达2的测量噪声的标准差，I2为2阶单位矩阵；
步骤1-3，通过雷达1获得目标初始时刻的状态，包含目标的位置和速度信息，得到第0个周期的目标状态强度的高斯成分集合为
Figure FDA00023419811900000212
其中，
Figure FDA00023419811900000213
为0时刻第i个高斯成分的权重，
Figure FDA00023419811900000214
为0时刻第i个高斯成分的状态向量，
Figure FDA00023419811900000215
为0时刻第i个高斯成分的方差阵，J0为0时刻高斯成分的数量；
步骤1-1中，状态转移矩阵Fk和噪声方差阵Qk分别如下：
Figure FDA00023419811900000216
其中，T为测量周期，
Figure FDA00023419811900000217
为系统噪声的标准差，F、Q为中间变量；
步骤2中，得到k时刻雷达1的预测高斯成分集合为
Figure FDA00023419811900000218
其中，
Figure FDA00023419811900000219
为k时刻第i个预测高斯成分的权重，
Figure FDA00023419811900000220
为k时刻第i个预测高斯成分的状态向量，
Figure FDA00023419811900000221
为k时刻第i个预测高斯成分的方差阵，Jk|k-1为k时刻预测高斯成分的数量；
步骤3包括如下步骤：
步骤3-1，根据雷达1的预测高斯成分集合计算预测测量值
Figure FDA00023419811900000222
预测测量误差协方差阵
Figure FDA0002341981190000031
增益矩阵
Figure FDA0002341981190000032
和估计误差协方差阵
Figure FDA0002341981190000033
公式如下：
Figure FDA0002341981190000034
其中，I为单位矩阵；
步骤3-2，对k时刻雷达1的每一个测量值
Figure FDA0002341981190000035
对预测高斯成分进行测量更新得到多目标状态集合矩
Figure FDA0002341981190000036
表示为：
Figure FDA0002341981190000037
其中，
Figure FDA0002341981190000038
表示将第i个预测高斯成分作为以x为变量、以
Figure FDA0002341981190000039
为均值、以
Figure FDA00023419811900000310
为方差阵的高斯分布，
Figure FDA00023419811900000311
为高斯分布的权重；
Figure FDA00023419811900000312
表示用第j个测量值更新第i个预测高斯成分后得到的更新后的以x为变量、以
Figure FDA00023419811900000313
为均值、以
Figure FDA00023419811900000314
为方差阵的高斯分布，其权重为
Figure FDA00023419811900000315
Figure FDA00023419811900000316
其中，
Figure FDA00023419811900000317
为第i个高斯成分的状态，
Figure FDA00023419811900000318
为k时刻第l个预测高斯成分的权重，
Figure FDA00023419811900000319
为第i个高斯成分的方差阵，kk(zj)为第j个测量值的杂波强度函数，
Figure FDA00023419811900000320
为雷达1由第1至k测量周期所有测量构成的集合，
Figure FDA00023419811900000321
为雷达1检测到目标的概率；
步骤3-3，得到测量更新后的高斯成分集合
Figure FDA00023419811900000322
其中，
Figure FDA00023419811900000323
为k时刻第i个更新后高斯成分的权重，
Figure FDA00023419811900000324
为k时刻第i个更新后高斯成分的状态向量，
Figure FDA00023419811900000325
为k时刻第i个更新后高斯成分的方差阵，Jk|k为k时刻更新后高斯成分的数量；
步骤4包括如下步骤：
步骤4-1，删除测量更新后的高斯成分集合
Figure FDA0002341981190000041
中权重充分小的高斯项，即删除权重
Figure FDA0002341981190000042
的高斯项，其中
Figure FDA0002341981190000043
为裁减门限；
步骤4-2，将距离充分小的高斯成分合并成一个，即将
Figure FDA0002341981190000044
的高斯项合并成一个，其中
Figure FDA0002341981190000045
为合并门限，距离di,j定义为
Figure FDA0002341981190000046
高斯成分合并方法如下：
Figure FDA0002341981190000047
其中L为所有可合并高斯成分的指标集合，
Figure FDA0002341981190000048
为k时刻第i个合并后高斯成分的权重，
Figure FDA0002341981190000049
为k时刻第i个合并后高斯成分的状态向量，
Figure FDA00023419811900000410
为k时刻第i个合并后高斯成分的方差阵，j为k时刻L中可合并高斯成分的下标；
步骤4-3，得到雷达1修剪融合后的高斯成分集合为
Figure FDA00023419811900000411
其中，
Figure FDA00023419811900000412
为修剪融合后的第i个高斯成分，
Figure FDA00023419811900000413
为其状态，
Figure FDA00023419811900000414
为其协方差阵，Jk,1为雷达1修剪融合后高斯成分的个数，为进行序贯操作，将该高斯成分集合作为雷达2的预测高斯成分集合；
步骤5包括如下步骤：
步骤5-1，根据雷达2的预测高斯成分集合计算预测测量值
Figure FDA00023419811900000415
预测测量误差协方差阵
Figure FDA00023419811900000416
增益矩阵
Figure FDA00023419811900000417
估计误差协方差阵
Figure FDA00023419811900000418
公式如下：
Figure FDA00023419811900000419
步骤5-2，对k时刻雷达2的每一个测量值
Figure FDA00023419811900000420
对预测高斯成分进行测量更新得到多目标状态集合矩
Figure FDA0002341981190000051
表示为：
Figure FDA0002341981190000052
其中，
Figure FDA0002341981190000053
其中，kk(zj)为杂波强度函数，
Figure FDA0002341981190000054
为雷达2由第1至k测量周期所有测量构成的集合，
Figure FDA0002341981190000055
为雷达2检测到目标的概率；
步骤5-3，得到测量更新后的高斯成分集合
Figure FDA0002341981190000056
步骤6包括如下步骤：
步骤6-1，删除测量更新后的高斯成分集合
Figure FDA0002341981190000057
中权重充分小的高斯项，即删除权重
Figure FDA0002341981190000058
的高斯项，其中
Figure FDA0002341981190000059
为裁减门限；
步骤6-2，将距离充分小的高斯成分合并成一个，即将
Figure FDA00023419811900000510
的高斯成分合并成一个，其中
Figure FDA00023419811900000511
为合并门限，距离di,j定义为
Figure FDA00023419811900000512
高斯成分合并方法如下：
Figure FDA00023419811900000513
其中L为所有可合并高斯成分的指标集合；
步骤6-3，得到雷达2修剪融合后的高斯成分集合为
Figure FDA00023419811900000514
其中，
Figure FDA00023419811900000515
为修剪融合后的第i个高斯成分，
Figure FDA00023419811900000516
为其状态，
Figure FDA00023419811900000517
为其协方差阵，Jk,2为雷达2修剪融合后高斯成分的个数。
2.根据权利要求1所述的方法，其特征在于，步骤7包括如下步骤：
步骤7-1，对于步骤4得到的雷达1的修剪融合后的高斯成分集合I1，即
Figure FDA0002341981190000061
中的每一个元素，从步骤6得到的雷达2修剪融合后的高斯成分集合I2，即
Figure FDA0002341981190000062
中找出与其距离最近的元素，设
Figure FDA0002341981190000063
为最近距离，M为维持融合门限，若dj＞M，则将I1中元素
Figure FDA0002341981190000064
添加到I2中，否则，通过两个元素中的权值项进行比较，将两个元素融合为一个元素，融合的具体过程为：
若I1中元素权值不小于0.5而I2中元素权值小于0.5，则用I1中元素替换I2中元素，若I1中元素权值小于0.5而I2中元素权值不小于0.5，则不进行操作，若两元素权值均小于或均不小于0.5，则进行加权平均，公式如下：
Figure FDA0002341981190000065
其中
Figure FDA0002341981190000066
为I1中元素，
Figure FDA0002341981190000067
为I2中元素，w1，w2和m为中间变量；
步骤7-2，得到维持融合后的高斯成分集合为
Figure FDA0002341981190000068
其中，
Figure FDA0002341981190000069
为维持融合后第i个高斯成分的权重，
Figure FDA00023419811900000610
为其状态，
Figure FDA00023419811900000611
为其协方差阵，Jk为维持融合后高斯成分的数量。
CN201810203353.XA 2018-03-13 2018-03-13 一种双雷达修正序贯高斯混合概率假设密度滤波方法 Active CN108732564B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810203353.XA CN108732564B (zh) 2018-03-13 2018-03-13 一种双雷达修正序贯高斯混合概率假设密度滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810203353.XA CN108732564B (zh) 2018-03-13 2018-03-13 一种双雷达修正序贯高斯混合概率假设密度滤波方法

Publications (2)

Publication Number Publication Date
CN108732564A CN108732564A (zh) 2018-11-02
CN108732564B true CN108732564B (zh) 2020-04-17

Family

ID=63940370

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810203353.XA Active CN108732564B (zh) 2018-03-13 2018-03-13 一种双雷达修正序贯高斯混合概率假设密度滤波方法

Country Status (1)

Country Link
CN (1) CN108732564B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109525220B (zh) * 2018-12-10 2019-08-30 中国人民解放军国防科技大学 具有航迹关联和提取能力的高斯混合cphd滤波方法
CN110376582B (zh) * 2019-01-24 2022-10-04 西安电子科技大学 自适应gm-phd的机动目标跟踪方法
CN111736145B (zh) * 2020-06-28 2022-04-19 电子科技大学 一种基于高斯混合概率假设密度滤波的多机动目标多普勒雷达跟踪方法
CN117724087B (zh) * 2024-02-07 2024-05-28 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103324835B (zh) * 2013-05-30 2016-09-28 深圳大学 概率假设密度滤波器目标信息的保持方法及信息保持系统
CN106022340A (zh) * 2016-05-17 2016-10-12 南京理工大学 一种改进的高斯混合势概率假设密度滤波方法
CN106896352A (zh) * 2017-04-17 2017-06-27 电子科技大学 一种基于随机集理论的多雷达异步数据分布式融合方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103324835B (zh) * 2013-05-30 2016-09-28 深圳大学 概率假设密度滤波器目标信息的保持方法及信息保持系统
CN106022340A (zh) * 2016-05-17 2016-10-12 南京理工大学 一种改进的高斯混合势概率假设密度滤波方法
CN106896352A (zh) * 2017-04-17 2017-06-27 电子科技大学 一种基于随机集理论的多雷达异步数据分布式融合方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于GMPHD的雷达组网检测跟踪算法研究;赵温波等;《系统仿真学报》;20161130;第28卷(第11期);2804-2812 *
组网无源雷达变数目多目标跟踪算法;时银水等;《西安电子科技大学学报》;20100430;第37卷(第02期);218-223 *

Also Published As

Publication number Publication date
CN108732564A (zh) 2018-11-02

Similar Documents

Publication Publication Date Title
CN108732564B (zh) 一种双雷达修正序贯高斯混合概率假设密度滤波方法
CN109633589A (zh) 目标跟踪中基于多模型优化多假设的多目标数据关联方法
CN109886305B (zh) 一种基于gm-phd滤波的多传感器非顺序量测异步融合方法
CN108344981B (zh) 面向杂波的多传感器异步检测tsbf多目标跟踪方法
CN110542885A (zh) 一种复杂交通环境下的毫米波雷达目标跟踪方法
CN109298725B (zh) 一种基于phd滤波的群体机器人分布式多目标跟踪方法
CN116128932B (zh) 一种多目标跟踪方法
CN111829505A (zh) 多传感器航迹质量外推航迹融合方法
CN106934324A (zh) 基于简化多假设算法的雷达数据关联方法
CN114002667B (zh) 一种基于随机矩阵法的多近邻扩展目标跟踪算法
CN110596693A (zh) 一种迭代更新的多传感器gmphd自适应融合方法
CN110109095A (zh) 目标特征辅助多源数据的关联方法
CN111127523A (zh) 基于量测迭代更新的多传感器gmphd自适应融合方法
CN110376581B (zh) 基于高斯混合概率假设密度滤波器的显式多目标跟踪方法
Chong et al. Ground target tracking-a historical perspective
CN114202025B (zh) 一种多传感器数据融合方法
CN115061139A (zh) 一种面向智能驾驶车辆的多传感器融合方法及系统
CN109214432B (zh) 一种多传感器多目标联合检测、跟踪与分类方法
CN111274529B (zh) 一种鲁棒的高斯逆威沙特phd多扩展目标跟踪算法
CN104637070A (zh) 基于概率假设密度的目标数变化的视频跟踪算法
CN109190647B (zh) 一种有源无源数据融合方法
Ma et al. Radiation intensity Gaussian mixture PHD filter for close target tracking
CN111340853B (zh) 基于ospa迭代的多传感器gmphd自适应融合方法
Blasch Modeling Intent for a target tracking and identification Scenario
CN115220002B (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