Nothing Special   »   [go: up one dir, main page]

CN111722214A - 雷达多目标跟踪phd实现方法 - Google Patents

雷达多目标跟踪phd实现方法 Download PDF

Info

Publication number
CN111722214A
CN111722214A CN202010493946.1A CN202010493946A CN111722214A CN 111722214 A CN111722214 A CN 111722214A CN 202010493946 A CN202010493946 A CN 202010493946A CN 111722214 A CN111722214 A CN 111722214A
Authority
CN
China
Prior art keywords
target
radar
gaussian
tracking
moment
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.)
Granted
Application number
CN202010493946.1A
Other languages
English (en)
Other versions
CN111722214B (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.)
Kunming University of Science and Technology
Original Assignee
Kunming University of Science and 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 Kunming University of Science and Technology filed Critical Kunming University of Science and Technology
Priority to CN202010493946.1A priority Critical patent/CN111722214B/zh
Publication of CN111722214A publication Critical patent/CN111722214A/zh
Application granted granted Critical
Publication of CN111722214B publication Critical patent/CN111722214B/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/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
    • G01S13/723—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
    • G01S13/726—Multiple target tracking
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

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

本发明公开了一种雷达多目标跟踪PHD实现方法，具体包括以下步骤：S1，基于多个跟踪目标在k‑1时刻的状态，计算各目标k时刻的状态预测估计；S2，将k时刻的雷达量测数据转换为笛卡尔坐标系下的目标状态似然数据；S3，将一步预测的高斯分量与目标状态似然数据进行乘积混合滤波；S4，计算雷达跟踪区域内所有目标的高斯分量和最终权重；S5，舍弃最终权重小于枝剪阈值的高斯分量，合并分布距离小于合并阈值的高斯分量；S6，令跟踪目标数目取所有高斯分量的权重和，目标状态取权重大于0.5的高斯分量的均值，迭代进行下一轮滤波直至跟踪结束；本发明的数据融合过程简单，抗杂波性能优异，对跟踪目标的状态估计精度高，目标数目估计准确。

Description

雷达多目标跟踪PHD实现方法
技术领域
本发明属于雷达多目标跟踪技术领域，特别是涉及一种雷达多目标跟踪PHD实现方法。
背景技术
随着传感技术、信息融合理论的飞速发展，雷达跟踪在军事和民用方面都得到了广泛的应用，例如空中监视、导弹防御、无人驾驶和气象监控等，然而随着跟踪环境的复杂化，特别是进入21世纪以后，隐身、电子对抗等新技术的日益完善，常规的单目标雷达跟踪方法已不足以处理这些问题。
上世纪五十年代，多目标跟踪技术(Multiple target tracking，MTT)被首次提出，最近六十年里，MTT在雷达跟踪领域取得了很多显著的成果，雷达跟踪多个目标时，主要问题在于目标存在新生、衍生和消亡等现象，雷达经常受到临近物体、山体环境、气象和电磁造成的影响，产生虚假量测；二十一世纪初国外学者提出通过随机有限集(Randomfinite set，RFS)描述多目标状态集和传感器观测集，利用最优贝叶斯滤波技术实现关联不确定、检测不确定和虚警杂波背景下的时变数目目标跟踪，由于多目标贝叶斯递推过程中包含多重积分，一般情况下采用多目标后验概率密度的一阶矩即概率假设密度(Probability Hypothesis Density，PHD)来代替多目标后验概率密度来进行计算。
与传统的雷达多目标跟踪方法相比，PHD滤波方法可以简单有效的描述虚警杂波、雷达跟踪目标的出现和消失等情况，不需要对雷达数据和目标状态进行复杂的关联计算，降低了运算量，具有优异的抗杂波性能，实现了复杂环境下数目未知的雷达多目标跟踪。
然而，PHD作为新兴的滤波技术，在实际应用中存在诸多不足，如实现困难、稳定性较差、不能直接处理非线性系统等；在跟踪监视多个目标时，雷达会受到杂波干扰产生虚假量测，其量测数据为雷达与目标的直线距离、雷达与目标的倾斜角度和俯仰角度，量测数据与笛卡尔坐标系下的目标位置构成一组非线性函数，现在多是将非线性滤波算法与PHD结合，这样处理计算量较大，且目标数目估计精度低、目标跟踪轨迹不稳定，特别是跟踪区域内目标数目越多时，雷达跟踪效果越差。
发明内容
本发明的目的在于提供一种雷达多目标跟踪PHD实现方法，以完成对复杂环境下雷达监测的多个目标的有效滤波跟踪，其处理过程简明高效、目标数目估计精度高、目标跟踪轨迹稳定、跟踪效果较好。
本发明所采用的技术方案是，雷达多目标跟踪PHD实现方法，具体包括以下步骤：
S1、已知k-1时刻雷达跟踪区域内多个目标的初始位置和运动模型，使用随机变量固定点采样非线性变换高斯近似方法，对k-1时刻雷达跟踪区域内各目标的每个高斯分量进行状态预测估计，得到雷达跟踪区域内各目标在k时刻的一步预测估计PHDDk|k-1(x)：
Dk|k-1(x)＝DS,k|k-1(x)+γk|k-1(x)
其中DS,k|k-1(x)为从k-1时刻到k时刻仍然存活的目标S的概率假设密度，γk|k-1(x)为从k-1时刻到k时刻新生的目标γ的概率假设密度；
将仍然存活的目标的高斯分量和新生目标的高斯分量统一为k时刻一步预测高斯分量，即：
Figure BDA0002522123530000021
其中Jk|k-1为基于k-1时刻各目标状态的高斯分量一步预测得到的k时刻高斯分量总数，i为表示k时刻一步预测高斯分量数目的变量，
Figure BDA0002522123530000022
为k时刻第i个高斯分量的权重，
Figure BDA0002522123530000023
表示目标状态x服从高斯分布，
Figure BDA0002522123530000024
分别为k时刻第i个高斯分量的均值和协方差；
S2、以雷达为坐标原点建立笛卡尔坐标系，使用Unscented变换将k时刻的雷达量测数据转换为笛卡尔坐标系下高斯形态的目标状态似然数据，得到雷达量测值zkj目标状态似然数据的均值
Figure BDA0002522123530000025
和协方差
Figure BDA0002522123530000026
S3、对目标状态似然数据与一步预测得到的高斯分量进行乘积融合滤波，得到融合后各高斯分量的似然权重、均值和协方差；
S4、使用PHD滤波器计算k时刻雷达漏检目标的高斯分量，得到k时刻雷达跟踪区域内所有目标的高斯分量Dk|k(x)和最终权重
Figure BDA0002522123530000027
S5、设置枝剪阈值和合并阈值，舍弃最终权重
Figure BDA0002522123530000031
小于枝剪阈值的高斯分量，合并分布距离db,e小于合并阈值的高斯分量；
S6、经过枝剪和合并后，令跟踪目标数目取所有高斯分量的权重和，目标状态取权重大于0.5的高斯分量的均值，迭代进行下一轮滤波直至跟踪结束。
进一步的，所述S1中使用的随机变量固定点采样非线性变换高斯近似方法，为Unscented变换、Cabuture变换或高斯厄米特数值积分。
进一步的，所述S2使用Unscented变换进行数据转换的过程如下：
以k时刻的第j个雷达量测值zkj为均值，以雷达量测噪声为协方差，按照Unscented变换规则生成sigma采样点
Figure BDA0002522123530000032
和采样点权值
Figure BDA0002522123530000033
利用反函数计算采样点并求出目标状态似然数据的均值和协方差：
Figure BDA0002522123530000034
Figure BDA0002522123530000035
Figure BDA0002522123530000036
其中n为雷达量测数据的维度，
Figure BDA0002522123530000037
为采样点
Figure BDA0002522123530000038
计算得到的样本点，h-1为雷达量测函数的反函数，
Figure BDA0002522123530000039
为雷达量测值zkj目标状态似然数据的均值，
Figure BDA00025221235300000310
为第a个样本点用于计算均值时的权重，
Figure BDA00025221235300000311
为雷达量测值zkj目标状态似然数据的协方差，
Figure BDA00025221235300000312
为第a个样本点用于计算协方差时的权重，T为转置；
对雷达量测值状态似然数据进行扩维处理：
Figure BDA00025221235300000313
Figure BDA00025221235300000314
Figure BDA00025221235300000315
为扩维后的
Figure BDA00025221235300000316
Figure BDA00025221235300000317
为扩维后的
Figure BDA00025221235300000318
F1 T为扩维矩阵F1的转置。
进一步的，所述S3中乘积融合滤波过程是：
Figure BDA00025221235300000319
Figure BDA0002522123530000041
Figure BDA0002522123530000042
其中
Figure BDA0002522123530000043
分别为第i个高斯分量与zkj的目标状态似然数据融合得到的高斯分量的似然权重、均值和协方差，
Figure BDA0002522123530000044
为
Figure BDA0002522123530000045
的反函数，
Figure BDA0002522123530000046
为
Figure BDA0002522123530000047
的反函数。
进一步的，所述S4中高斯分量Dk|k(x)和最终权重
Figure BDA0002522123530000048
的计算公式如下：
Figure BDA0002522123530000049
Figure BDA00025221235300000410
其中PD,k为k时刻雷达对目标的检测概率，PD,k取0.95～0.99，Zk为k时刻雷达量测值的集合，
Figure BDA00025221235300000411
为最终权重，即第i个高斯分量和第j个雷达量测值乘积融合得到的高斯分量在生成Dk|k(x)时的权重，λ为杂波泊松分布的均值，c(zkj)为k时刻第j个雷达量测值来自杂波的概率密度。
进一步的，所述S5中设置的枝剪阈值为0.0001，合并阈值为5，合并距离
Figure BDA00025221235300000412
k时刻雷达跟踪区域内所有目标的高斯分量的合并过程为：
Figure BDA00025221235300000413
Figure BDA00025221235300000414
Figure BDA00025221235300000415
其中b、e为表示k时刻雷达跟踪区域内所有目标的高斯分量数目的变量，b≠e，
Figure BDA00025221235300000416
为第b个高斯分量的协方差，
Figure BDA00025221235300000417
为
Figure BDA00025221235300000418
的反函数，wb,e为第b个高斯分量的权重
Figure BDA00025221235300000419
与第e个高斯分量的权重
Figure BDA00025221235300000420
合并后得到的权重，mb,e为第b个高斯分量的均值
Figure BDA00025221235300000421
与第e个高斯分量的均值
Figure BDA00025221235300000422
合并后得到的均值，pb,e为第b个高斯分量的协方差
Figure BDA00025221235300000423
与第e个高斯分量的协方差
Figure BDA0002522123530000051
合并后得到的协方差。
本发明的有益效果是：与现有非线性算法结合PHD的雷达跟踪滤波方法比较，本发明给出了一种全新的雷达多目标跟踪PHD实现方法，该方法直接对多目标先验信息和雷达量测信息进行乘积融合，无需将多目标先验信息转换为雷达量测数据，简化了数据融合过程；并通过计算残差进行滤波，抑制了预测过程噪声的影响，简明易于理解，抗杂波性能优异，特别是在雷达量测数据误差比系统误差小的情况下，状态估计精度高，目标数目估计准确。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例的流程简图。
图2是本发明实施例一次仿真的量测数据和真实轨迹图。
图3是本发明实施例一次仿真的GM-CK-PHD算法滤波轨迹和真实轨迹的比较图。
图4是本发明实施例一次仿真的GM-UK-PHD算法滤波轨迹和真实轨迹的比较图。
图5是本发明实施例一次仿真的本发明所得轨迹和真实轨迹的比较图。
图6是三种多目标跟踪方法的跟踪目标数目估计和真实目标数目的对比图。
图7是使用本发明方法和GM-CK-PHD、GM-UK-PHD得到的100次蒙特卡洛仿真的OSPA距离比较图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
参照图1，雷达多目标跟踪PHD实现方法具体包括以下步骤：
S1，在雷达多目标跟踪非线性系统中，采用混合高斯形式的概率假设密度描述雷达跟踪监视区域内的目标，使用高斯分量表示雷达跟踪监视区域内目标可能存在的状态；已知雷达跟踪区域内多个目标在k-1时刻的初始位置和运动模型，使用随机变量固定点采样非线性变换高斯近似方法，对k-1时刻雷达跟踪区域内各目标的每个高斯分量进行状态预测估计，得到雷达跟踪区域内各目标k时刻的状态预测估计PHDDk|k-1(x)：
Dk|k-1(x)＝DS,k|k-1(x)+γk|k-1(x)
随机变量固定点采样非线性变换高斯近似方法为Unscented变换、Cabuture变换或高斯厄米特数值积分，DS,k|k-1(x)为从k-1时刻到k时刻仍然存活的目标S的概率假设密度，γk|k-1(x)为从k-1时刻到k时刻新生的目标γ的概率假设密度；
设PS,k为k时刻目标存活的概率，N(x；m,p)为目标状态x服从高斯分布，m和p分别为目标状态的均值和协方差，w为目标状态的权重，即目标存在的强度，则：
Figure BDA0002522123530000061
Figure BDA0002522123530000062
其中l为表示k-1时刻存活目标的高斯分量数目变量，JS,k-1为k-1时刻存活目标高斯分量总数，
Figure BDA0002522123530000063
和
Figure BDA0002522123530000064
分别为k-1时刻存活目标第l个高斯分量的权重、均值和协方差，q为表示新生目标高斯分量数目变量，Jγ,k为k时刻新生目标高斯分量总数，
Figure BDA0002522123530000065
和
Figure BDA0002522123530000066
分别为k时刻新生目标第q个高斯分量的权重、均值和协方差；
为了方便表达，将仍然存活的目标的高斯分量和新生目标的高斯分量统一为k时刻一步预测高斯分量，即：
Figure BDA0002522123530000067
其中Jk|k-1为基于k-1时刻各目标状态的高斯分量一步预测得到的k时刻高斯分量总数，i为表示k时刻一步预测高斯分量的数目变量，
Figure BDA0002522123530000068
分别为k时刻第i个高斯分量的权重、均值和协方差；
S2，由于雷达量测数据为极坐标数据，被跟踪目标状态数据为笛卡尔坐标数据，为实现乘积融合，以雷达为坐标原点建立笛卡尔坐标系，使用Unscented变换将k时刻的雷达量测数据转换为笛卡尔坐标系下高斯形态的状态似然数据，得到雷达量测值zkj目标状态似然数据的均值
Figure BDA0002522123530000069
和协方差
Figure BDA00025221235300000610
所述使用Unscented变换进行数据转换的过程如下：
以k时刻的第j个雷达量测值zkj为均值，以雷达量测噪声为协方差，按照Unscented变换规则生成sigma采样点
Figure BDA0002522123530000071
和采样点权值
Figure BDA0002522123530000072
即对雷达量测数据进行固定规则采样；
利用反函数计算样本点并求出雷达量测数据状态似然的一阶矩、二阶矩，即均值和协方差：
Figure BDA0002522123530000073
Figure BDA0002522123530000074
Figure BDA0002522123530000075
其中n为雷达量测数据的维度，
Figure BDA0002522123530000076
为采样点
Figure BDA0002522123530000077
计算得到的样本点，h-1为雷达量测函数的反函数，在二维平面，h-1为笛卡尔坐标系下被跟踪目标的位置坐标，其横坐标与纵坐标分别为雷达与被跟踪目标的直线距离乘以角度数据的正弦值和余弦值，
Figure BDA0002522123530000078
为雷达量测值zkj目标状态似然数据的均值，
Figure BDA0002522123530000079
为第a个样本点用于计算均值时的权重，
Figure BDA00025221235300000710
为雷达量测值zkj目标状态似然数据的协方差，
Figure BDA00025221235300000711
为第a个样本点用于计算协方差时的权重，T为转置；
由于雷达量测值反向估计得到的目标状态似然数据只有位置信息，即状态向量的维度数为2，与一步预测估计数据维度不匹配，所以需要对反向估计得到的目标状态似然数据进行扩维，即在无信息维度处置零，扩维过程如下：
Figure BDA00025221235300000712
Figure BDA00025221235300000713
Figure BDA00025221235300000714
为扩维后的
Figure BDA00025221235300000715
Figure BDA00025221235300000716
为扩维后的
Figure BDA00025221235300000717
F1 T为扩维矩阵F1的转置；
S3，根据概率似然乘积规则，对雷达量测数据得出的目标状态似然数据与一步预测得到的多目标状态高斯分量进行乘积融合滤波，得到融合后多目标状态高斯分量的权重、均值和协方差，所述乘积融合滤波过程如下：
Figure BDA00025221235300000718
Figure BDA0002522123530000081
Figure BDA0002522123530000082
其中
Figure BDA0002522123530000083
分别为第i个高斯分量与zkj的目标状态似然数据融合得到的高斯分量的似然权重、均值和协方差，
Figure BDA0002522123530000084
为
Figure BDA0002522123530000085
的反函数，
Figure BDA0002522123530000086
为
Figure BDA0002522123530000087
的反函数；
S4，通过PHD滤波器计算k时刻雷达漏检目标高斯分量，得到k时刻雷达跟踪区域内所有目标的高斯分量Dk|k(x)和最终权重
Figure BDA0002522123530000088
以消除杂波干扰、突出跟踪目标的权重，具体公式为：
Figure BDA0002522123530000089
Figure BDA00025221235300000810
其中PD,k为k时刻雷达对目标的检测概率，一般取0.95至0.99之间，Zk为k时刻的雷达量测值的集合，
Figure BDA00025221235300000811
为最终权重，即第i个高斯分量和第j个雷达量测值乘积融合得到的高斯分量在生成Dk|k(x)时的权重，设定杂波强度服从泊松分布，杂波空间分布的概率密度为监视区域的均匀分布，λ为泊松分布的均值，c(zkj)为k时刻第j个雷达量测值来自杂波的概率密度；
S5，为防止迭代造成的高斯分量过多，计算量增大，设定枝剪阈值为0.0001，合并阈值为5，舍弃最终权重
Figure BDA00025221235300000812
小于枝剪阈值的高斯分量，对分布距离db,e小于合并阈值的高斯分量进行合并，
Figure BDA00025221235300000813
k时刻雷达跟踪区域内所有目标的高斯分量的合并过程为：
Figure BDA00025221235300000814
Figure BDA00025221235300000815
Figure BDA00025221235300000816
其中b、e为表示k时刻雷达跟踪区域内所有目标的高斯分量数目的变量，b≠e，wb,e为第b个高斯分量的权重
Figure BDA0002522123530000091
与第e个高斯分量的权重
Figure BDA0002522123530000092
合并后得到的权重，mb,e为第b个高斯分量的均值
Figure BDA0002522123530000093
与第e个高斯分量的均值
Figure BDA0002522123530000094
合并后得到的均值，pb,e为第b个高斯分量的协方差
Figure BDA0002522123530000095
与第e个高斯分量的协方差
Figure BDA0002522123530000096
合并后得到的协方差；
S6，经过枝剪和合并后，目标数目估计取所有高斯分量的权重和，目标状态取权重大于0.5的高斯分量的均值，迭代进行下一轮滤波更新，直至跟踪结束。
实施例
使用二维雷达对多目标进行跟踪，并使用本发明对多目标跟踪过程进行状态预测，本实施例中目标做近似匀速转弯运动，其运动方程和雷达量测方程如下：
xk＝Fxk-1+Gwk
Figure BDA0002522123530000097
其中xk为k时刻目标的状态变量，状态变量
Figure BDA0002522123530000098
Figure BDA0002522123530000099
采样周期t＝1s，
Figure BDA00025221235300000910
为过程噪声，
Figure BDA00025221235300000911
为雷达量测距离误差，
Figure BDA00025221235300000912
雷达量测角度误差，取σw＝0.2m/s，σr＝0.06m，σθ＝0.1rad，角速度ω＝0.4rad/s；以雷达为坐标原点建立笛卡尔坐标系，xk和
Figure BDA00025221235300000913
分别为X轴方向的位置和速度，yk和
Figure BDA00025221235300000914
分别为Y轴方向的位置和速度，Zk为k时刻的雷达量测数据集，θk为k时刻的角度量测值，rk为k时刻的距离量测值。
参照图1，使用本发明跟踪多目标的具体步骤为：
假设k-1时刻PHD的混合高斯参数为
Figure BDA00025221235300000915
其中JS,k-1为k-1时刻存活目标的高斯分量总数，
Figure BDA00025221235300000916
和
Figure BDA00025221235300000917
分别为k-1时刻第l个高斯分量的权重、均值和协方差，结合多目标先验信息和Unscented非线性变换方法求出k时刻的一步预测PHD：
Dk|k-1(x)＝DS,k|k-1(x)+γk|k-1(x)
其中DS,k|k-1(x)为从k-1时刻到k时刻仍然存活的目标S的概率假设密度，γk|k-1(x)为k时刻新生目标γ的概率假设密度，
Figure BDA0002522123530000101
多目标跟踪中目标存活概率PS,k＝0.99，
Figure BDA0002522123530000102
获得k时刻的雷达量测数据集合Zk＝{zk1,zk2,…,zkM(k)}，其中M(k)为k时刻总的观测数据，采用Unscented变换方法生成样本点和权值：
Figure BDA0002522123530000103
Figure BDA0002522123530000104
按照Unscented变换规则计算出状态似然一阶矩、二阶矩：
Figure BDA0002522123530000105
Figure BDA0002522123530000106
Figure BDA0002522123530000107
由于雷达观测向量反向估计得到的目标状态似然只有位置信息，即状态向量的维度数为2，与多目标状态一步预测估计数据维度不匹配，所以对反向估计得到的目标状态似然数据进行扩维处理，扩维矩阵为：
Figure BDA0002522123530000108
扩维处理如下：
Figure BDA0002522123530000109
Figure BDA00025221235300001010
乘积融合：使用概率似然乘积融合规则将一步预测高斯分量与扩维后的目标状态似然数据进行乘积融合：
Figure BDA00025221235300001011
Figure BDA00025221235300001012
Figure BDA0002522123530000111
Figure BDA0002522123530000112
为第i个高斯分量与zkj的目标状态似然数据融合得到的高斯分量的似然权重，权重代表目标强度即目标存在的可能性。
现有方法中权值都是通过将目标状态信息转换成观测信息后，结合真实观测数据计算似然函数得到，而本发明不需要对目标状态信息进行转换，而是将观测信息转换为目标状态信息，减少了过程噪声的影响，特别是在观测噪声较小的情况下，本发明计算得到的权值更准确。
在PHD滤波器框架下更新k时刻的PHD，得到k时刻雷达跟踪区域内所有目标的高斯分量Dk|k(x)，
Figure BDA0002522123530000113
PD,k为k时刻雷达对目标的检测概率，PD,k＝0.99，
Figure BDA0002522123530000114
杂波在观测空间上服从均匀分布c(z)，数目服从参数λ＝5的泊松分布。
最后设置剪枝阈值C＝1e-4，合并阈值为5，对k时刻雷达跟踪区域内所有目标的高斯分量进行枝剪、合并，提取目标状态和数目，储存后验估计信息，进入下一步迭代。
高斯分量最多为Jmax＝100，共采取100个时刻，目标真实运动参数如表1所示，其中目标3为目标2在50s时的衍生目标，新生目标的强度都为0.1，仿真选取OSPA距离作为算法性能的评价指标，OSPA距离越小，算法精度越高，数目估计越准确，其中参数p＝2，c＝8。
表1目标真实运动参数
目标 起始时刻/s 消失时刻/s 目标初始状态
1 1 60 (-65m,1.5m/s,70m,-1.5m/s)
2 22 80 (-50m,1.5m/s,60m,-1.5m/s)
3 50 80 目标2在50s时的衍生
4 41 100 (65m,-1.5m/s,-55m,-1.5m/s)
5 51 100 (-65m,1.6m/s,-45m,-1.6m/s)
为了验证本发明的有效性，仿真实例将本发明(GM-PP-PHD)与容积卡尔曼(CKF)实现的GM-CK-PHD、无迹卡尔曼(UKF)实现的GM-UK-PHD做了对比，图2为量测数据和目标真实轨迹图，图3、图4和图5为三种算法的跟踪轨迹图，由图3、图4可以看出GM-CK-PHD和GM-UK-PHD虽然能跟踪目标的真实轨迹，但是对目标位置进行估计时数据有部分缺失，且在目标移动轨迹以外获得了数个干扰数据，而由图5可知本发明能跟踪目标的真实轨迹，且目标位置数据较为完整、无干扰数据，即本发明的跟踪轨迹稳定，跟踪效果较好；图6为目标数目估计比较图，由图6可知GM-CK-PHD和GM-UK-PHD跟踪的目标数目距离真实目标数目有较大的起伏，而本发明跟踪的目标数目与真实目标个数较为接近、起伏较少，本发明对目标数目的跟踪估计精度较高。
图7为100次蒙特卡洛仿真不同算法的OSPA距离比较图，表2为三种算法100次蒙特卡洛仿真的平均时间。
表2三种算法的平均耗时
Figure BDA0002522123530000121
由图7可以看出，GM-PP-PHD滤波器的OSPA距离比另外两种算法都小，且一直比较稳定，即本发明提出的GM-PP-PHD实现方法目标位置和数目估计更准确，且从表2可知，相对于经典算法本发明并未增加计算量，具有较高的工程应用价值。
本说明书中的各个实施例均采用相关的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

Claims (6)

1.雷达多目标跟踪PHD实现方法，其特征在于，具体包括以下步骤：
S1、已知k-1时刻雷达跟踪区域内多个目标的初始位置和运动模型，使用随机变量固定点采样非线性变换高斯近似方法，对k-1时刻雷达跟踪区域内各目标的每个高斯分量进行状态预测估计，得到雷达跟踪区域内各目标在k时刻的一步预测估计PHDDk|k-1(x)：
Dk|k-1(x)＝DS,k|k-1(x)+γk|k-1(x)
其中DS,k|k-1(x)为从k-1时刻到k时刻仍然存活的目标S的概率假设密度，γk|k-1(x)为从k-1时刻到k时刻新生的目标γ的概率假设密度；
将仍然存活的目标的高斯分量和新生目标的高斯分量统一为k时刻一步预测高斯分量，即：
Figure FDA0002522123520000011
其中Jk|k-1为基于k-1时刻各目标状态的高斯分量一步预测得到的k时刻高斯分量总数，i为表示k时刻一步预测高斯分量数目的变量，
Figure FDA0002522123520000012
为k时刻第i个高斯分量的权重，
Figure FDA0002522123520000013
表示目标状态x服从高斯分布，
Figure FDA0002522123520000014
分别为k时刻第i个高斯分量的均值和协方差；
S2、以雷达为坐标原点建立笛卡尔坐标系，使用Unscented变换将k时刻的雷达量测数据转换为笛卡尔坐标系下高斯形态的目标状态似然数据，得到雷达量测值zkj目标状态似然数据的均值
Figure FDA0002522123520000015
和协方差
Figure FDA0002522123520000016
S3、对目标状态似然数据与一步预测得到的高斯分量进行乘积融合滤波，得到融合后各高斯分量的似然权重、均值和协方差；
S4、使用PHD滤波器计算k时刻雷达漏检目标的高斯分量，得到k时刻雷达跟踪区域内所有目标的高斯分量Dk|k(x)和最终权重
Figure FDA0002522123520000017
S5、设置枝剪阈值和合并阈值，舍弃最终权重
Figure FDA0002522123520000018
小于枝剪阈值的高斯分量，合并分布距离db,e小于合并阈值的高斯分量；
S6、经过枝剪和合并后，令跟踪目标数目取所有高斯分量的权重和，目标状态取权重大于0.5的高斯分量的均值，迭代进行下一轮滤波直至跟踪结束。
2.根据权利要求1所述的雷达多目标跟踪PHD实现方法，其特征在于，所述S1中使用的随机变量固定点采样非线性变换高斯近似方法，为Unscented变换、Cabuture变换或高斯厄米特数值积分。
3.根据权利要求1所述的雷达多目标跟踪PHD实现方法，其特征在于，所述S2使用Unscented变换进行数据转换的过程如下：
以k时刻的第j个雷达量测值zkj为均值，以雷达量测噪声为协方差，按照Unscented变换规则生成sigma采样点
Figure FDA0002522123520000021
和采样点权值
Figure FDA0002522123520000022
利用反函数计算采样点并求出目标状态似然数据的均值和协方差：
Figure FDA0002522123520000023
Figure FDA0002522123520000024
Figure FDA0002522123520000025
其中n为雷达量测数据的维度，
Figure FDA0002522123520000026
为采样点
Figure FDA0002522123520000027
计算得到的样本点，h-1为雷达量测函数的反函数，
Figure FDA0002522123520000028
为雷达量测值zkj目标状态似然数据的均值，
Figure FDA0002522123520000029
为第a个样本点用于计算均值时的权重，
Figure FDA00025221235200000210
为雷达量测值zkj目标状态似然数据的协方差，
Figure FDA00025221235200000211
为第a个样本点用于计算协方差时的权重，T为转置；
对雷达量测值状态似然数据进行扩维处理：
Figure FDA00025221235200000212
Figure FDA00025221235200000213
Figure FDA00025221235200000214
为扩维后的
Figure FDA00025221235200000215
Figure FDA00025221235200000216
为扩维后的
Figure FDA00025221235200000217
F1 T为扩维矩阵F1的转置。
4.根据权利要求1所述的雷达多目标跟踪PHD实现方法，其特征在于，所述S3中乘积融合滤波过程是：
Figure FDA00025221235200000218
Figure FDA00025221235200000219
Figure FDA0002522123520000031
其中
Figure FDA0002522123520000032
分别为第i个高斯分量与zkj的目标状态似然数据融合得到的高斯分量的似然权重、均值和协方差，
Figure FDA0002522123520000033
为
Figure FDA0002522123520000034
的反函数，
Figure FDA0002522123520000035
为
Figure FDA0002522123520000036
的反函数。
5.根据权利要求1所述的雷达多目标跟踪PHD实现方法，其特征在于，所述S4中高斯分量Dk|k(x)和最终权重
Figure FDA0002522123520000037
的计算公式如下：
Figure FDA0002522123520000038
Figure FDA0002522123520000039
其中PD,k为k时刻雷达对目标的检测概率，PD,k取0.95～0.99，Zk为k时刻雷达量测值的集合，
Figure FDA00025221235200000310
为最终权重，即第i个高斯分量和第j个雷达量测值乘积融合得到的高斯分量在生成Dk|k(x)时的权重，λ为杂波泊松分布的均值，c(zkj)为k时刻第j个雷达量测值来自杂波的概率密度。
6.根据权利要求1所述的雷达多目标跟踪PHD实现方法，其特征在于，所述S5中设置的枝剪阈值为0.0001，合并阈值为5，合并距离
Figure FDA00025221235200000311
k时刻雷达跟踪区域内所有目标的高斯分量的合并过程为：
Figure FDA00025221235200000312
Figure FDA00025221235200000313
Figure FDA00025221235200000314
其中b、e为表示k时刻雷达跟踪区域内所有目标的高斯分量数目的变量，b≠e，
Figure FDA00025221235200000315
为第b个高斯分量的协方差，
Figure FDA00025221235200000316
为
Figure FDA00025221235200000317
的反函数，wb,e为第b个高斯分量的权重
Figure FDA00025221235200000318
与第e个高斯分量的权重
Figure FDA00025221235200000319
合并后得到的权重，mb,e为第b个高斯分量的均值
Figure FDA00025221235200000320
与第e个高斯分量的均值
Figure FDA00025221235200000321
合并后得到的均值，pb,e为第b个高斯分量的协方差
Figure FDA00025221235200000322
与第e个高斯分量的协方差
Figure FDA0002522123520000041
合并后得到的协方差。
CN202010493946.1A 2020-06-03 2020-06-03 雷达多目标跟踪phd实现方法 Active CN111722214B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010493946.1A CN111722214B (zh) 2020-06-03 2020-06-03 雷达多目标跟踪phd实现方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010493946.1A CN111722214B (zh) 2020-06-03 2020-06-03 雷达多目标跟踪phd实现方法

Publications (2)

Publication Number Publication Date
CN111722214A true CN111722214A (zh) 2020-09-29
CN111722214B CN111722214B (zh) 2024-01-30

Family

ID=72565649

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010493946.1A Active CN111722214B (zh) 2020-06-03 2020-06-03 雷达多目标跟踪phd实现方法

Country Status (1)

Country Link
CN (1) CN111722214B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112911249A (zh) * 2021-03-23 2021-06-04 浙江大华技术股份有限公司 目标对象的跟踪方法、装置、存储介质及电子装置
CN113093097A (zh) * 2021-03-18 2021-07-09 南京航空航天大学 一种使用互质阵列的概率假设密度doa跟踪的方法
CN113219452A (zh) * 2021-05-07 2021-08-06 电子科技大学 未知视域下的分布式多雷达联合配准与多目标跟踪方法
CN113238218A (zh) * 2021-04-22 2021-08-10 中山大学 基于phd滤波的临近空间高超声速目标跟踪方法
CN113537299A (zh) * 2021-06-23 2021-10-22 西北工业大学 基于期望最大化高斯混合缩减的分布式贝叶斯滤波器目标跟踪方法
CN114739383A (zh) * 2022-03-04 2022-07-12 中国人民解放军空军工程大学 一种基于多目标不确定性改进的gm-phd滤波方法
CN116047495A (zh) * 2023-03-31 2023-05-02 昆明理工大学 一种用于三坐标雷达的状态变换融合滤波跟踪方法
CN117724087A (zh) * 2024-02-07 2024-03-19 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106526585A (zh) * 2016-10-26 2017-03-22 中国人民解放军空军工程大学 基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法
CN109284677A (zh) * 2018-08-16 2019-01-29 昆明理工大学 一种贝叶斯滤波目标跟踪算法
CN110320512A (zh) * 2019-07-09 2019-10-11 大连海事大学 一种基于带标签的gm-phd平滑滤波多目标跟踪方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106526585A (zh) * 2016-10-26 2017-03-22 中国人民解放军空军工程大学 基于高斯粒子势概率假设密度滤波的目标检测前跟踪方法
CN109284677A (zh) * 2018-08-16 2019-01-29 昆明理工大学 一种贝叶斯滤波目标跟踪算法
CN110320512A (zh) * 2019-07-09 2019-10-11 大连海事大学 一种基于带标签的gm-phd平滑滤波多目标跟踪方法

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113093097A (zh) * 2021-03-18 2021-07-09 南京航空航天大学 一种使用互质阵列的概率假设密度doa跟踪的方法
CN112911249A (zh) * 2021-03-23 2021-06-04 浙江大华技术股份有限公司 目标对象的跟踪方法、装置、存储介质及电子装置
CN112911249B (zh) * 2021-03-23 2023-03-24 浙江大华技术股份有限公司 目标对象的跟踪方法、装置、存储介质及电子装置
CN113238218B (zh) * 2021-04-22 2023-10-20 中山大学 基于phd滤波的临近空间高超声速目标跟踪方法
CN113238218A (zh) * 2021-04-22 2021-08-10 中山大学 基于phd滤波的临近空间高超声速目标跟踪方法
CN113219452A (zh) * 2021-05-07 2021-08-06 电子科技大学 未知视域下的分布式多雷达联合配准与多目标跟踪方法
CN113219452B (zh) * 2021-05-07 2022-05-31 电子科技大学 未知视域下的分布式多雷达联合配准与多目标跟踪方法
CN113537299A (zh) * 2021-06-23 2021-10-22 西北工业大学 基于期望最大化高斯混合缩减的分布式贝叶斯滤波器目标跟踪方法
CN114739383A (zh) * 2022-03-04 2022-07-12 中国人民解放军空军工程大学 一种基于多目标不确定性改进的gm-phd滤波方法
CN114739383B (zh) * 2022-03-04 2024-05-14 中国人民解放军空军工程大学 一种基于多目标不确定性改进的gm-phd滤波方法
CN116047495B (zh) * 2023-03-31 2023-06-02 昆明理工大学 一种用于三坐标雷达的状态变换融合滤波跟踪方法
CN116047495A (zh) * 2023-03-31 2023-05-02 昆明理工大学 一种用于三坐标雷达的状态变换融合滤波跟踪方法
CN117724087A (zh) * 2024-02-07 2024-03-19 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法
CN117724087B (zh) * 2024-02-07 2024-05-28 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Also Published As

Publication number Publication date
CN111722214B (zh) 2024-01-30

Similar Documents

Publication Publication Date Title
CN111722214A (zh) 雷达多目标跟踪phd实现方法
CN109002835B (zh) 一种基于最大熵模糊聚类的粒子滤波数据关联方法
CN110503071B (zh) 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法
CN101975575B (zh) 基于粒子滤波的被动传感器多目标跟踪方法
Yang et al. Linear-time joint probabilistic data association for multiple extended object tracking
CN111711432B (zh) 一种基于ukf和pf混合滤波的目标跟踪算法
CN111127523A (zh) 基于量测迭代更新的多传感器gmphd自适应融合方法
CN111708013B (zh) 一种距离坐标系目标跟踪滤波方法
CN113673565B (zh) 多传感器gm-phd自适应序贯融合多目标跟踪方法
CN115204212A (zh) 一种基于stm-pmbm滤波算法的多目标跟踪方法
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
CN109115228A (zh) 一种基于加权最小二乘容积卡尔曼滤波的目标定位方法
CN108320302B (zh) 基于随机超曲面的CBMeMBer多目标跟踪方法
Yun et al. Variational Bayesian based adaptive PDA filter in scenarios with unknown detection probability and heavy-tailed process noise
CN111736144B (zh) 一种仅用距离观测的机动转弯目标状态估计方法
CN116047495B (zh) 一种用于三坐标雷达的状态变换融合滤波跟踪方法
CN112379350A (zh) 智能车辆毫米波雷达多目标跟踪方法、装置及设备
Chen et al. IMM tracking of a 3D maneuvering target with passive TDOA system
CN116224320B (zh) 一种极坐标系下处理多普勒量测的雷达目标跟踪方法
CN111340853A (zh) 基于ospa迭代的多传感器gmphd自适应融合方法
Wang et al. Infrared/laser multi-sensor fusion and tracking based on the multi-scale model
Lai et al. Adaptive multinoulli-based Kalman filter with randomly unknown delayed and lost measurements
CN115327503A (zh) 基于高斯粒子滤波的固定单站无源定位方法和相关装置
Fu et al. Observation bootstrapping CDK-GMPHD filter based on consensus fusion strategy
CN108981707B (zh) 基于时差量测箱粒子phd的被动跟踪多目标方法

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