CN107462882B - 一种适用于闪烁噪声的多机动目标跟踪方法及系统 - Google Patents

一种适用于闪烁噪声的多机动目标跟踪方法及系统 Download PDF

Info

Publication number
CN107462882B
CN107462882B CN201710804689.7A CN201710804689A CN107462882B CN 107462882 B CN107462882 B CN 107462882B CN 201710804689 A CN201710804689 A CN 201710804689A CN 107462882 B CN107462882 B CN 107462882B
Authority
CN
China
Prior art keywords
target
distribution
gamma
current moment
gaussian
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.)
Expired - Fee Related
Application number
CN201710804689.7A
Other languages
English (en)
Other versions
CN107462882A (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.)
Shenzhen University
Original Assignee
Shenzhen 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 Shenzhen University filed Critical Shenzhen University
Priority to CN201710804689.7A priority Critical patent/CN107462882B/zh
Publication of CN107462882A publication Critical patent/CN107462882A/zh
Priority to PCT/CN2018/071580 priority patent/WO2019047455A1/zh
Application granted granted Critical
Publication of CN107462882B publication Critical patent/CN107462882B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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
    • G01S11/00Systems for determining distance or velocity not using reflection or reradiation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/25Fusion techniques
    • G06F18/251Fusion techniques of input or preprocessed data

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Measurement Of Radiation (AREA)

Abstract

本发明适用于多传感器信息融合领域,提供了一种适用于闪烁噪声的多机动目标跟踪方法与系统,通过预测步骤、更新步骤、融合步骤、生成步骤以及提取步骤,用t分布建模闪烁噪声,应用变分贝叶斯方法近似地求出不同模型下观测噪声的协方差、目标状态以及自由度的联合概率密度,有效解决了闪烁噪声下多机动目标的跟踪问题,提高了多机动目标的跟踪精度。

Description

一种适用于闪烁噪声的多机动目标跟踪方法及系统
技术领域
本发明属于多传感器信息融合技术领域,尤其涉及一种适用于闪烁噪声的多机动目标跟踪方法及系统。
背景技术
闪烁噪声下的概率假设密度滤波器是适用于闪烁噪声多目标跟踪的有效方法,该滤波器能够对闪烁噪声下的非机动目标进行有效的跟踪,但该滤波器不能对在不同模型之间转换的机动目标进行跟踪。
发明内容
本发明提供了一种适用于闪烁噪声的多机动目标跟踪方法及系统,旨在解决闪烁噪声下多机动目标的有效跟踪问题。
为解决上述技术问题,本发明是这样实现的,本发明提供了一种适用于闪烁噪声的多机动目标跟踪方法,所述方法包括:
步骤1、利用前一时刻各个目标的运动模型、以及高斯伽玛混合分布和存在概率产生伽玛分布的形状参数和尺度参数;根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、模型间的马尔可夫转移概率、以及伽玛分布的形状参数和尺度参数,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率;
步骤2、根据当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率,利用变分贝叶斯方法对当前时刻的测量数据进行序贯处理,得到当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率;
步骤3、对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率分别进行融合,得到当前时刻已存在目标的运动模型、高斯伽玛混合分布和存在概率;
步骤4、利用当前时刻的测量数据生成新生目标的高斯伽玛混合分布,为所述新生目标指定存在概率,并选取匀速运动模型作为新生目标的运动模型,将所述新生目标的运动模型、高斯伽玛混合分布及存在概率分别与所述当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,得到当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率;
步骤5、从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标,并将裁减后余下的目标的运动模型、高斯伽玛混合分布和存在概率作为滤波器下一次递归的输入;
从所述裁减后余下的目标中提取存在概率大于第二阈值的目标,并将提取出的目标的高斯伽玛混合分布作为当前时刻的输出,将输出的高斯伽玛混合分布中,高斯分布的均值作为当前时刻目标的状态估计。
进一步地,所述步骤1包括:
以k-1表示前一时刻,k表示当前时刻,tk-1表示前一时刻的时间,tk表示当前时刻的时间,ri,k-1表示前一时刻目标i的运动模型标识,且1≤ri,k-1≤Mr,ri,k表示当前时刻目标i的模型标识,且ri,k=1,2,…,Mr,其中,Mr表示运动模型的总数目;当前时刻的测量噪声服从ξ维的学生t分布,以S(yj,k;Hkxi,k,Ri,ki,k)表示当前时刻测量的概率密度函数,其中,yj,k表示当前时刻的一个测量,Hk表示测量矩阵,xi,k表示目标的状态,Ri,k表示精度矩阵,νi,k表示t分布的自由度,且
Figure BDA0001402419060000021
为精度矩阵Ri,k主对角线上的各个元素;
已知k-1时刻目标i的运动模型为ri,k-1,则在该模型ri,k-1下的高斯伽玛混合分布为
Figure BDA0001402419060000031
存在概率为ρi,k-1(ri,k-1);其中,N表示高斯分布,g表示伽玛分布,xi,k-1表示k-1时刻目标i的状态向量,mi,k-1(ri,k-1)和Pi,k-1(ri,k-1)分别表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中高斯分布的均值和协方差,
Figure BDA0001402419060000032
表示k-1时刻精度矩阵
Figure BDA0001402419060000033
的对角线元素,νi,k-1表示k-1时刻t分布的自由度,
Figure BDA0001402419060000034
和γi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的形状参数,
Figure BDA0001402419060000035
Figure BDA00014024190600000310
表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的尺度参数,ξ表示状态向量的维数,i=1,2,...,Nk-1,Nk-1表示前一时刻目标的总数;
根据k-1时刻目标i的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、以及模型间的马尔可夫转移概率,得到当前时刻目标i在模型ri,k下的预测高斯伽玛混合分布为
Figure BDA0001402419060000036
预测存在概率为ρi,k|k-1(ri,k);其中,mi,k|k-1(ri,k)=Fk-1(ri,k)mi,k-1(ri,k-1),
Figure BDA0001402419060000037
Figure BDA0001402419060000038
γi,k|k-1(ri,k)=ργγi,k-1(ri,k-1),ηi,k|k-1(ri,k)=ρηηi,k-1(ri,k-1),ρi,k|k-1(ri,k)=ps,k(tk-tk-1)tk|k-1(ri,k|ri,k-1i,k-1(ri,k-1),Fk-1(ri,k)表示k-1时刻目标i在模型ri,k下的状态转移矩阵,Qk-1(ri,k)表示k-1时刻目标i在模型ri,k下的过程噪声协方差矩阵,上标T表示矩阵或向量的转置,ρα,ρβ,ργ,ρη表示传递因子,取值范围为(0,1],ps,k(tk-tk-1)表示目标的幸存概率,且
Figure BDA0001402419060000039
tk|k-1(ri,k|ri,k-1)表示模型间的马尔可夫转移概率,T表示采样周期,δ为预定的常数,l=1,...,ξ,ξ表示状态向量的维数。
进一步地,设当前时刻的测量集为yk=(y1,k,...,yM,k),其中,M为当前时刻的测量数据的总数,则所述步骤2包括:
以当前时刻各个目标在Mr个运动模型下的预测高斯伽玛混合分布和预测存在概率,作为当前时刻各个目标在相应运动模型下的初始高斯伽玛混合分布和初始存在概率,则目标i在模型ri,k下的初始高斯伽玛混合分布取为
Figure BDA0001402419060000041
初始存在概率取为
Figure BDA0001402419060000042
其中i=1,2,...,Nk-1,ri,k=1,2,…,Mr
Figure BDA0001402419060000043
Figure BDA0001402419060000044
利用变分贝叶斯方法对第1个至第M个当前时刻的测量数据进行序贯处理,具体包括以下:
设第j个测量数据处理前,目标i在模型ri,k下的高斯伽玛混合分布及存在概率分别为
Figure BDA0001402419060000045
Figure BDA0001402419060000046
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr,j=1,2,…,M;
根据变分贝叶斯方法,计算得到第j个测量更新时,目标i在模型ri,k下的高斯伽玛 混合分布为 存在概率分别为其中,
Figure BDA0001402419060000051
;表示伽玛函数,tr表示矩阵的迹, 表示均值向量,表示协方差矩阵,表示滤波器增益,伽玛分布 的形状参数为和伽玛分布的尺度参数为和精度矩阵为其中,Hk为模型ri,k的观测矩阵,PD,k为目标的检测概 率,λc,k为杂波密度,yj,k为当前时刻接收到的第j个测量数据,上标T表示矩阵或向量的转 置,ξ为状态向量的总维数;
Figure BDA0001402419060000061
则第j个测量数据处理后,目标i在模型ri,k下的的高斯伽玛混合分布为
Figure BDA0001402419060000062
目标i的存在概率为
Figure BDA0001402419060000063
其中
Figure BDA0001402419060000064
Figure BDA0001402419060000065
Figure BDA0001402419060000066
则第j个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布为
Figure BDA0001402419060000067
目标i的存在概率为
Figure BDA0001402419060000068
其中,
Figure BDA0001402419060000069
Figure BDA00014024190600000610
在第M个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布和存在概率分别为
Figure BDA00014024190600000611
Figure BDA00014024190600000613
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr
将第M个测量数据处理后,各个目标在不同模型下的高斯伽玛混合分布及存在概率分别作为当前时刻各个目标的更新高斯伽玛混合分布和存在概率,则所述当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布为
Figure BDA00014024190600000612
更新存在概率为
Figure BDA0001402419060000071
其中i=1,…,Nk-1,ri,k=1,2,…,Mr
Figure BDA0001402419060000072
Figure BDA0001402419060000073
Figure BDA0001402419060000074
进一步地,所述步骤3包括:
对当前时刻目标i在Mr个运动模型下的更新的高斯伽玛混合分布和存在概率分别进行融合,得到当前时刻已存在目标q的运动模型为rq,k、高斯伽玛混合分布为
Figure BDA0001402419060000075
和存在概率为ρq,k(rq,k),其中,q为当前时刻已存在目标的索引号,且q=1,…,Nk-1
Figure BDA0001402419060000076
Figure BDA0001402419060000077
伽玛分布的形状参数为
Figure BDA0001402419060000078
Figure BDA0001402419060000079
伽玛分布的尺度参数为
Figure BDA00014024190600000710
Figure BDA00014024190600000711
进一步地,所述步骤4包括:
选取匀速运动模型作为新生目标的运动模型,设此时新生目标的运动模型为rj,k=1,其中,j=1,…,M;
利用M个测量数据生成当前时刻新生目标的高斯伽玛混合分布为
Figure BDA00014024190600000712
其中,j=1,…,M,
Figure BDA0001402419060000081
为第j个新生目标的协方差,
Figure BDA0001402419060000082
由第j个测量数据yj,k=[xj,k yj,k]T产生,且
Figure BDA00014024190600000811
Figure BDA0001402419060000084
为第j个新生目标分布函数中伽玛分布的形状参数,
Figure BDA0001402419060000085
Figure BDA0001402419060000086
为第j个新生目标分布函数中伽玛分布的尺度参数;
指定当前时刻各新生目标的存在概率为
Figure BDA0001402419060000087
其中,j=1,…,M,ρυ为所指定的存在概率;
将所述当前时刻已存在目标的运动模型与所述当前时刻新生目标的运动模型合并,得到当前时刻各个目标的运动模型为
Figure BDA0001402419060000088
其中,Nk=Nk-1+M;
将所述当前时刻已存在目标的高斯伽玛混合分布与所述当前时刻新生目标的高斯伽玛混合分布进行合并,得到当前时刻各个目标的高斯伽玛混合分布为
Figure BDA0001402419060000089
将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并,得到所述当前时刻各个目标的存在概率为
Figure BDA00014024190600000810
为解决上述技术问题,本发明还提供了一种适用于闪烁噪声的多机动目标跟踪系统,所述系统包括:
预测模块,用于利用前一时刻各个目标的运动模型、以及高斯伽玛混合分布和存在概率产生伽玛分布的形状参数和尺度参数;根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、模型间的马尔可夫转移概率、以及伽玛分布的形状参数和尺度参数,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率;
更新模块,用于根据当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率,利用变分贝叶斯方法对当前时刻的测量数据进行序贯处理,得到当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率;
融合模块,用于对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率分别进行融合,得到当前时刻已存在目标的运动模型、高斯伽玛混合分布和存在概率;
生成模块,用于利用当前时刻的测量数据生成新生目标的高斯伽玛混合分布,为所述新生目标指定存在概率,并选取匀速运动模型作为新生目标的运动模型,将所述新生目标的运动模型、高斯伽玛混合分布及存在概率分别与所述当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,得到当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率;
提取模块,用于从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标,并将裁减后余下的目标的运动模型、高斯伽玛混合分布和存在概率作为滤波器下一次递归的输入;从所述裁减后余下的目标中提取存在概率大于第二阈值的目标,并将提取出的目标的高斯伽玛混合分布作为当前时刻的输出,将输出的高斯伽玛混合分布中,高斯分布的均值作为当前时刻目标的状态估计。
本发明与现有技术相比,有益效果在于:
本发明所提供的方法通过用t分布建模闪烁噪声,应用变分贝叶斯方法近似地求出不同模型下观测噪声的协方差、目标状态以及自由度的联合概率密度,有效解决了闪烁噪声下多机动目标的跟踪问题,提高了多机动目标的跟踪精度。
附图说明
图1是本发明第一个实施例提供的一种适用于闪烁噪声的多机动目标跟踪方法示意图;
图2是本发明实施例提供的传感器经过70个扫描周期得到的测量数据示意图;
图3是在闪烁噪声下采用JMS-PHD滤波器的多目标跟踪方法处理得到的结果示意图;
图4是在闪烁噪声下采用本发明所提供的多机动目标跟踪方法处理得到的结果示意图;
图5是本发明实施例提供的按照本发明所提供的多目标跟踪方法与JMS-PHD滤波方法在经过100次实验得到的平均OSPA距离对比示意图;
图6是本发明第二个实施例提供一种适用于闪烁噪声的多目标跟踪系统的结构示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
作为本发明的第一个实施例,如图1所示,本发明提供的一种适用于闪烁噪声的多机动目标跟踪方法,该方法包括以下步骤:
步骤S101:在马尔科夫链控制的不同系统模型中,利用前一时刻各个目标的运动模型、以及高斯伽玛混合分布和存在概率产生伽玛分布的形状参数和尺度参数;根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、模型间的马尔可夫转移概率、以及伽玛分布的形状参数和尺度参数,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率。步骤S101具体包括以下步骤:
以k-1表示前一时刻,k表示当前时刻,tk-1表示前一时刻的时间,tk表示当前时刻的时间,ri,k-1表示前一时刻目标i的运动模型标识,且1≤ri,k-1≤Mr,ri,k表示当前时刻目标i的模型标识,且ri,k=1,2,…,Mr,其中,Mr表示运动模型的总数目;当前时刻的测量噪声服从ξ维的学生t分布,以S(yj,k;Hkxi,k,Ri,ki,k)表示当前时刻测量的概率密度函数,其中,yj,k表示当前时刻的一个测量,Hk表示测量矩阵,xi,k表示目标的状态,Ri,k表示精度矩阵,νi,k表示t分布的自由度,且
Figure BDA0001402419060000111
为精度矩阵Ri,k主对角线上的各个元素;
已知k-1时刻目标i的运动模型为ri,k-1,则在该模型ri,k-1下的高斯伽玛混合分布为
Figure BDA0001402419060000112
存在概率为ρi,k-1(ri,k-1);其中,N表示高斯分布,g表示伽玛分布,xi,k-1表示k-1时刻目标i的状态向量,mi,k-1(ri,k-1)和Pi,k-1(ri,k-1)分别表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中高斯分布的均值和协方差,
Figure BDA0001402419060000113
表示k-1时刻精度矩阵
Figure BDA0001402419060000114
的对角线元素,νi,k-1表示k-1时刻t分布的自由度,
Figure BDA0001402419060000115
和γi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的形状参数,
Figure BDA0001402419060000116
和ηi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的尺度参数,ξ表示状态向量的维数,i=1,2,...,Nk-1,Nk-1表示前一时刻目标的总数;
根据k-1时刻目标i的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、以及模型间的马尔可夫转移概率,得到当前时刻目标i在模型ri,k下的预测高斯伽玛混合分布为
Figure BDA0001402419060000117
预测存在概率为ρi,k|k-1(ri,k);其中,mi,k|k-1(ri,k)=Fk-1(ri,k)mi,k-1(ri,k-1),
Figure BDA0001402419060000118
Figure BDA0001402419060000119
γi,k|k-1(ri,k)=ργγi,k-1(ri,k-1),ηi,k|k-1(ri,k)=ρηηi,k-1(ri,k-1),ρi,k|k-1(ri,k)=ps,k(tk-tk-1)tk|k-1(ri,k|ri,k-1i,k-1(ri,k-1),Fk-1(ri,k)表示k-1时刻目标i在模型ri,k下的状态转移矩阵,Qk-1(ri,k)表示k-1时刻目标i在模型ri,k下的过程噪声协方差矩阵,上标T表示矩阵或向量的转置,ρα,ρβ,ργ,ρη表示传递因子,取值范围为(0,1],ps,k(tk-tk-1)表示目标的幸存概率,且
Figure BDA0001402419060000121
tk|k-1(ri,k|ri,k-1)表示模型间的马尔可夫转移概率,T表示采样周期,δ为预定的常数,l=1,...,ξ,ξ表示状态向量的维数。
步骤S102:根据当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率,利用变分贝叶斯方法对当前时刻的测量数据进行序贯处理,得到当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率。设当前时刻的测量集为yk=(y1,k,...,yM,k),其中,M为当前时刻的测量数据的总数,则此时步骤S102具体包括以下步骤:
以当前时刻各个目标在Mr个运动模型下的预测高斯伽玛混合分布和预测存在概率,作为当前时刻各个目标在相应运动模型下的初始高斯伽玛混合分布和初始存在概率,则目标i在模型ri,k下的初始高斯伽玛混合分布取为
Figure BDA0001402419060000122
初始存在概率取为
Figure BDA0001402419060000123
其中i=1,2,...,Nk-1,ri,k=1,2,…,Mr
Figure BDA0001402419060000124
Figure BDA0001402419060000125
利用变分贝叶斯方法对第1个至第M个当前时刻的测量数据进行序贯处理,具体包括以下:
设第j个测量数据处理前,目标i在模型ri,k下的高斯伽玛混合分布及存在概率分别为
Figure BDA0001402419060000131
Figure BDA0001402419060000132
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr,j=1,2,…,M;
根据变分贝叶斯方法,计算得到第j个测量更新时,目标i在模型ri,k下的高斯伽玛 混合分布为 存在概率分别为其中,
Figure BDA0001402419060000135
;表示伽玛函数,tr表示矩阵的迹, 表示均值向量,表示协方差矩阵,表示滤波器增益,伽玛分布 的形状参数为和伽玛分布的尺度参数为和精度矩阵为其中,Hk为模型ri,k的观测矩阵,PD,k为目标的检测概 率,λc,k为杂波密度,yj,k为当前时刻接收到的第j个测量数据,上标T表示矩阵或向量的转 置,ξ为状态向量的总维数;
Figure BDA0001402419060000146
则第j个测量数据处理后,目标i在模型ri,k下的的高斯伽玛混合分布为
Figure BDA0001402419060000147
目标i的存在概率为
Figure BDA0001402419060000148
其中
Figure BDA0001402419060000149
Figure BDA00014024190600001410
Figure BDA00014024190600001411
则第j个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布为
Figure BDA00014024190600001412
目标i的存在概率为
Figure BDA00014024190600001413
其中,
Figure BDA00014024190600001414
Figure BDA00014024190600001415
在第M个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布和存在概率分别为
Figure BDA0001402419060000151
Figure BDA0001402419060000159
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr
将第M个测量数据处理后,各个目标在不同模型下的高斯伽玛混合分布及存在概率分别作为当前时刻各个目标的更新高斯伽玛混合分布和存在概率,则所述当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布为
Figure BDA0001402419060000152
更新存在概率为
Figure BDA0001402419060000153
其中i=1,…,Nk-1,ri,k=1,2,…,Mr
Figure BDA0001402419060000154
Figure BDA0001402419060000155
Figure BDA00014024190600001510
步骤S103:对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率分别进行融合,得到当前时刻已存在目标的运动模型、高斯伽玛混合分布和存在概率。步骤S103具体包括以下步骤:
对当前时刻目标i在Mr个运动模型下的更新的高斯伽玛混合分布和存在概率分别进行融合,得到当前时刻已存在目标q的运动模型为rq,k、高斯伽玛混合分布为
Figure BDA0001402419060000156
和存在概率为ρq,k(rq,k),其中,q为当前时刻已存在目标的索引号,且q=1,…,Nk-1
Figure BDA0001402419060000157
Figure BDA0001402419060000158
伽玛分布的形状参数为
Figure BDA0001402419060000161
Figure BDA0001402419060000162
伽玛分布的尺度参数为
Figure BDA0001402419060000163
Figure BDA0001402419060000164
步骤S104:利用当前时刻的测量数据生成新生目标的高斯伽玛混合分布,为所述新生目标指定存在概率,并选取匀速运动模型作为新生目标的运动模型,将所述新生目标的运动模型、高斯伽玛混合分布及存在概率分别与所述当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,得到当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率。步骤S104具体包括以下步骤:
选取匀速运动模型作为新生目标的运动模型,设此时新生目标的运动模型为rj,k=1,其中,j=1,…,M。该匀速运动模型为多种模型中的其中一种类型(可以定义为模型1)。
利用M个测量数据生成当前时刻新生目标的高斯伽玛混合分布为
Figure BDA0001402419060000165
其中,j=1,…,M,
Figure BDA0001402419060000166
为第j个新生目标的协方差,
Figure BDA0001402419060000167
由第j个测量数据yj,k=[xj,k yj,k]T产生,且
Figure BDA00014024190600001614
Figure BDA00014024190600001610
为第j个新生目标分布函数中伽玛分布的形状参数,
Figure BDA00014024190600001611
Figure BDA00014024190600001612
为第j个新生目标分布函数中伽玛分布的尺度参数;
指定当前时刻各新生目标的存在概率为
Figure BDA00014024190600001613
其中,j=1,…,M,ρυ为所指定的存在概率;
将所述当前时刻已存在目标的运动模型与所述当前时刻新生目标的运动模型合并,得到当前时刻各个目标的运动模型为
Figure BDA0001402419060000171
其中,Nk=Nk-1+M;
将所述当前时刻已存在目标的高斯伽玛混合分布与所述当前时刻新生目标的高斯伽玛混合分布进行合并,得到当前时刻各个目标的高斯伽玛混合分布为
Figure BDA0001402419060000172
将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并,得到所述当前时刻各个目标的存在概率为
Figure BDA0001402419060000173
步骤S105:从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标,并将裁减后余下的目标的运动模型、高斯伽玛混合分布和存在概率作为滤波器下一次递归的输入;从所述裁减后余下的目标中提取存在概率大于第二阈值的目标,并将提取出的目标的高斯伽玛混合分布作为当前时刻的输出,将输出的高斯伽玛混合分布中,高斯分布的均值作为当前时刻目标的状态估计。在本实施例中,经过多次的仿真实验得出,将第一阈值设定为10-3,将第二阈值设定为0.5,可以达到较好的效果。
下面结合图2至图5对本实施例进行进一步地解释:
在本实施例中,考虑二维空间[-400(m),400(m)]×[-400(m),400(m)]中运动的机动目标。目标的状态由位置和速度构成,表示为
Figure BDA0001402419060000174
其中x和y分别表示位置分量,
Figure BDA0001402419060000175
Figure BDA0001402419060000176
分别表示速度分量,上标T表示向量的转置,状态转移矩阵为
Figure BDA0001402419060000181
过程噪声方差矩阵为
Figure BDA0001402419060000182
Δtk=tk-tk-1为当前时刻与前一时刻的时间差,σv为过程噪声标准差;目标的运动模型包括ri,k=1,ri,k=2,ri,k=3,模型1是一个匀速运动模型,该模型1的状态转移矩阵和协方差矩阵分别为Fk-1(ri,k=1)=F(ω=0s-1)和Qk-1(ri,k=1)=Q(σv=1ms-2);模型2是一个匀转弯模型,该模型2的状态转移矩阵和协方差矩阵分别为Fk-1(ri,k=2)=F(ω=5°s-1)和Qk-1(ri,k=2)=Q(σv=3ms-2);模型3是一个匀转弯模型,该模型3的状态转移矩阵和协方差矩阵分别为Fk-1(ri,k=3)=F(ω=-5°s-1)和Qk-1(ri,k=3)=Q(σv=3ms-2);不同运动模型之间的马尔科夫转移矩阵为
Figure BDA0001402419060000183
观测矩阵
Figure BDA0001402419060000184
观测噪声假设服从rk=10的t分布。
为了产生仿真数据,取幸存概率pS,k=1.0、检测概率pD,k=0.9、杂波密度λc,k=4.6875×10-6m-2。一次实验中传感器经过70个扫描周期得到的仿真观测数据如图2所示。本实施例在闪烁噪声下,分别采用JMS-PHD(跳变马尔科夫系统模型概率假设密度)滤波器的多目标跟踪方法和本发明所提供的适用于闪烁噪声的多机动目标跟踪方法对图2中的仿真数据分别进行了处理,在采用JMS-PHD滤波器进行处理时,相关参数设置为pS,k=1.0、pD,k=0.9、λc,k=4.6875×10-6m-2、第一阈值为10-3、第二阈值为0.5,传递因子ρα=ρβ=ργ=ρη=0.98,伽玛分布形状参数初始值
Figure BDA0001402419060000191
伽玛分布尺度参数初始值
Figure BDA0001402419060000192
在闪烁噪声下,采用JMS-PHD滤波器时的新生目标的权重为wυ=0.1,采用本发明实施例所提供的方法时的新生目标的存在概率为ρυ=0.1、新生目标的协方差
Figure BDA0001402419060000193
图3和图4分别为JMS-PHD滤波器方法所产生的结果和本发明实施例提供的多机动目标跟踪方法所产生的结果,通过两图的对比结果可以直观的得出,采用本发明所提供的目标跟踪方法所得到的运动轨迹更加连续、完整、精确。通过本发明实施例所提供的方法与现有的闪烁噪声下的JMS-PHD滤波器对图2的仿真数据分别进行处理,经过100次Monte Carlo实验,得到平均OSPA(Optimal Subpattern Assignment,最优亚模式分配)距离,如图5所示。通过图中对比可以看出,本发明的多目标跟踪方法在闪烁噪声下对于关联不确定、检测不确定的多机动目标跟踪可获得更为精确和可靠的目标状态估计、其滤波后的OSPA距离比闪烁噪声下的JMS-PHD滤波器滤波后的OSPA距离要小。
综上所述,本发明第一个实施例所提供的方法,根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率,以及模型间的马尔可夫转移概率,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布及预测存在概率,根据各个目标在不同模型中预测的高斯伽玛混合分布和存在概率,利用变分贝叶斯方法序贯处理当前时刻接收到的各测量得到各目标在不同模型下的更新高斯伽玛混合分布和存在概率,对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和存在概率进行融合,得到当前时刻各个已存在目标的运动模型、高斯伽玛混合分布和存在概率,利用当前时刻的测量数据生成新目标的高斯伽玛混合分布,为所述新目标指定存在概率,并选取匀速运动模型作为新目标的运动模型,再将新生目标的运动模型、高斯伽玛混合分布及存在概率分别与当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,生成当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率,从而很好地解决了闪烁噪声下多机动目标的跟踪问题,可以达到有效跟踪的目的。
作为本发明的第二个实施例,如图6所示,本发明提供了一种适用于闪烁噪声的多机动目标跟踪系统,所述系统包括:
预测模块101,用于利用前一时刻各个目标的运动模型、以及高斯伽玛混合分布和存在概率产生伽玛分布的形状参数和尺度参数;根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、模型间的马尔可夫转移概率、以及伽玛分布的形状参数和尺度参数,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率。预测模块101具体用于:
以k-1表示前一时刻,k表示当前时刻,tk-1表示前一时刻的时间,tk表示当前时刻的时间,ri,k-1表示前一时刻目标i的运动模型标识,且1≤ri,k-1≤Mr,ri,k表示当前时刻目标i的模型标识,且ri,k=1,2,…,Mr,其中,Mr表示运动模型的总数目;当前时刻的测量噪声服从ξ维的学生t分布,以S(yj,k;Hkxi,k,Ri,ki,k)表示当前时刻测量的概率密度函数,其中,yj,k表示当前时刻的一个测量,Hk表示测量矩阵,xi,k表示目标的状态,Ri,k表示精度矩阵,νi,k表示t分布的自由度,且
Figure BDA0001402419060000201
为精度矩阵Ri,k主对角线上的各个元素;
已知k-1时刻目标i的运动模型为ri,k-1,则在该模型ri,k-1下的高斯伽玛混合分布为
Figure BDA0001402419060000202
存在概率为ρi,k-1(ri,k-1);其中,N表示高斯分布,g表示伽玛分布,xi,k-1表示k-1时刻目标i的状态向量,mi,k-1(ri,k-1)和Pi,k-1(ri,k-1)分别表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中高斯分布的均值和协方差,
Figure BDA0001402419060000203
表示k-1时刻精度矩阵
Figure BDA0001402419060000211
的对角线元素,νi,k-1表示k-1时刻t分布的自由度,
Figure BDA0001402419060000212
和γi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的形状参数,
Figure BDA0001402419060000213
Figure BDA0001402419060000214
表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的尺度参数,ξ表示状态向量的维数,i=1,2,...,Nk-1,Nk-1表示前一时刻目标的总数;
根据k-1时刻目标i的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、以及模型间的马尔可夫转移概率,得到当前时刻目标i在模型ri,k下的预测高斯伽玛混合分布为
Figure BDA0001402419060000215
预测存在概率为ρi,k|k-1(ri,k);其中,mi,k|k-1(ri,k)=Fk-1(ri,k)mi,k-1(ri,k-1),
Figure BDA0001402419060000216
Figure BDA0001402419060000217
γi,k|k-1(ri,k)=ργγi,k-1(ri,k-1),ηi,k|k-1(ri,k)=ρηηi,k-1(ri,k-1),ρi,k|k-1(ri,k)=ps,k(tk-tk-1)tk|k-1(ri,k|ri,k-1i,k-1(ri,k-1),Fk-1(ri,k)表示k-1时刻目标i在模型ri,k下的状态转移矩阵,Qk-1(ri,k)表示k-1时刻目标i在模型ri,k下的过程噪声协方差矩阵,上标T表示矩阵或向量的转置,ρα,ρβ,ργ,ρη表示传递因子,取值范围为(0,1],ps,k(tk-tk-1)表示目标的幸存概率,且
Figure BDA0001402419060000218
tk|k-1(ri,k|ri,k-1)表示模型间的马尔可夫转移概率,T表示采样周期,δ为预定的常数,l=1,...,ξ,ξ表示状态向量的维数。
更新模块102,用于根据当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率,利用变分贝叶斯方法对当前时刻的测量数据进行序贯处理,得到当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率。设当前时刻的测量集为yk=(y1,k,...,yM,k),其中,M为当前时刻测量的总数,则更新模块102具体用于:
以当前时刻各个目标在Mr个运动模型下的预测高斯伽玛混合分布和预测存在概率,作为当前时刻各个目标在相应运动模型下的初始高斯伽玛混合分布和初始存在概率,则目标i在模型ri,k下的初始高斯伽玛混合分布取为
Figure BDA0001402419060000221
初始存在概率取为
Figure BDA0001402419060000222
其中i=1,2,...,Nk-1,ri,k=1,2,…,Mr
Figure BDA0001402419060000223
Figure BDA0001402419060000228
利用变分贝叶斯方法对第1个至第M个当前时刻的测量数据进行序贯处理,具体包括以下:
设第j个测量数据处理前,目标i在模型ri,k下的高斯伽玛混合分布及存在概率分别为
Figure BDA0001402419060000224
Figure BDA0001402419060000225
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr,j=1,2,…,M;
根据变分贝叶斯方法,计算得到第j个测量更新时,目标i在模型ri,k下的高斯伽玛 混合分布为 存在概率分别为其中,
Figure BDA0001402419060000231
;表示伽玛函数,tr表示矩阵的迹, 表示均值向量,表示协方差矩阵,表示滤波器增益,伽玛分布 的形状参数为和伽玛分布的尺度参数为和精度矩阵为其中,Hk为模型ri,k的观测矩阵,PD,k为目标的检测概 率,λc,k为杂波密度,yj,k为当前时刻接收到的第j个测量数据,上标T表示矩阵或向量的转 置,ξ为状态向量的总维数;
Figure BDA0001402419060000241
则第j个测量数据处理后,目标i在模型ri,k下的的高斯伽玛混合分布为
Figure BDA0001402419060000242
目标i的存在概率为
Figure BDA0001402419060000243
其中
Figure BDA0001402419060000244
Figure BDA0001402419060000245
Figure BDA0001402419060000246
则第j个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布为
Figure BDA0001402419060000247
目标i的存在概率为
Figure BDA0001402419060000248
其中,
Figure BDA0001402419060000249
Figure BDA00014024190600002410
在第M个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布和存在概率分别为
Figure BDA00014024190600002411
Figure BDA00014024190600002412
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr
将第M个测量数据处理后,各个目标在不同模型下的高斯伽玛混合分布及存在概率分别作为当前时刻各个目标的更新高斯伽玛混合分布和存在概率,则所述当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布为
Figure BDA00014024190600002413
Figure BDA00014024190600002414
更新存在概率为
Figure BDA0001402419060000251
其中i=1,…,Nk-1,ri,k=1,2,…,Mr
Figure BDA0001402419060000252
Figure BDA0001402419060000253
Figure BDA0001402419060000254
融合模块103,用于对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率分别进行融合,得到当前时刻已存在目标的运动模型、高斯伽玛混合分布和存在概率。融合模块103具体用于:
对当前时刻目标i在Mr个运动模型下的更新的高斯伽玛混合分布和存在概率分别进行融合,得到当前时刻已存在目标q的运动模型为rq,k、高斯伽玛混合分布为
Figure BDA0001402419060000255
和存在概率为ρq,k(rq,k),其中,q为当前时刻已存在目标的索引号,且q=1,…,Nk-1
Figure BDA0001402419060000256
Figure BDA0001402419060000257
伽玛分布的形状参数为
Figure BDA0001402419060000258
Figure BDA0001402419060000259
伽玛分布的尺度参数为
Figure BDA00014024190600002510
Figure BDA00014024190600002511
生成模块104,用于利用当前时刻的测量数据生成新生目标的高斯伽玛混合分布,为所述新生目标指定存在概率,并选取匀速运动模型作为新生目标的运动模型,将所述新生目标的运动模型、高斯伽玛混合分布及存在概率分别与所述当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,得到当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率。生成模块104具体用于:
选取匀速运动模型作为新生目标的运动模型,设此时新生目标的运动模型为rj,k=1,其中,j=1,…,M;
利用M个测量数据生成当前时刻新生目标的高斯伽玛混合分布为
Figure BDA0001402419060000261
其中,j=1,…,M,
Figure BDA0001402419060000262
为第j个新生目标的协方差,
Figure BDA0001402419060000263
由第j个测量数据yj,k=[xj,k yj,k]T产生,且
Figure BDA00014024190600002612
Figure BDA0001402419060000266
为第j个新生目标分布函数中伽玛分布的形状参数,
Figure BDA0001402419060000267
Figure BDA0001402419060000268
为第j个新生目标分布函数中伽玛分布的尺度参数;
指定当前时刻各新生目标的存在概率为
Figure BDA0001402419060000269
其中,j=1,…,M,ρυ为所指定的存在概率;
将所述当前时刻已存在目标的运动模型与所述当前时刻新生目标的运动模型合并,得到当前时刻各个目标的运动模型为
Figure BDA00014024190600002610
其中,Nk=Nk-1+M;
将所述当前时刻已存在目标的高斯伽玛混合分布与所述当前时刻新生目标的高斯伽玛混合分布进行合并,得到当前时刻各个目标的高斯伽玛混合分布为
Figure BDA00014024190600002611
将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并,得到所述当前时刻各个目标的存在概率为
Figure BDA0001402419060000271
提取模块105,用于从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标,并将裁减后余下的目标的运动模型、高斯伽玛混合分布和存在概率作为滤波器下一次递归的输入;从所述裁减后余下的目标中提取存在概率大于第二阈值的目标,并将提取出的目标的高斯伽玛混合分布作为当前时刻的输出,将输出的高斯伽玛混合分布中,高斯分布的均值作为当前时刻目标的状态估计。
综上所述,本发明第二个实施例所提供的系统,通过预测模块、更新模块、融合模块、生成模块以及提取模块,用t分布建模闪烁噪声,应用变分贝叶斯方法近似地求出不同模型下观测噪声的协方差、目标状态以及自由度的联合概率密度,有效解决了闪烁噪声下多机动目标的跟踪问题,提高了多机动目标的跟踪精度。
以上所述仅为本发明的较佳实施例而已,并不用以限制发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。

Claims (8)

1.一种适用于闪烁噪声的多机动目标跟踪方法,其特征在于,所述方法包括:
步骤1、利用前一时刻各个目标的运动模型、以及高斯伽玛混合分布和存在概率产生伽玛分布的形状参数和尺度参数;根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、模型间的马尔可夫转移概率、以及伽玛分布的形状参数和尺度参数,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率;所述步骤1包括:
以k-1表示前一时刻,k表示当前时刻,tk-1表示前一时刻的时间,tk表示当前时刻的时间,ri,k-1表示前一时刻目标i的运动模型标识,且1≤ri,k-1≤Mr,ri,k表示当前时刻目标i的模型标识,且ri,k=1,2,…,Mr,其中,Mr表示运动模型的总数目;当前时刻的测量噪声服从ξ维的学生t分布,以S(yj,k;Hkxi,k,Ri,ki,k)表示当前时刻测量的概率密度函数,其中,yj,k表示当前时刻的一个测量,Hk表示测量矩阵,xi,k表示目标的状态,Ri,k表示精度矩阵,νi,k表示t分布的自由度,且
Figure FDA0002400948370000011
Figure FDA0002400948370000012
为精度矩阵Ri,k主对角线上的各个元素;
已知k-1时刻目标i的运动模型为ri,k-1,则在该模型ri,k-1下的高斯伽玛混合分布为
Figure FDA0002400948370000013
存在概率为ρi,k-1(ri,k-1);其中,N表示高斯分布,g表示伽玛分布,xi,k-1表示k-1时刻目标i的状态向量,mi,k-1(ri,k-1)和Pi,k-1(ri,k-1)分别表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中高斯分布的均值和协方差,
Figure FDA0002400948370000014
表示k-1时刻精度矩阵
Figure FDA0002400948370000015
的对角线元素,νi,k-1表示k-1时刻t分布的自由度,
Figure FDA0002400948370000016
和γi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的形状参数,
Figure FDA0002400948370000021
和ηi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的尺度参数,ξ表示状态向量的维数,i=1,2,...,Nk-1,Nk-1表示前一时刻目标的总数;
根据k-1时刻目标i的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、以及模型间的马尔可夫转移概率,得到当前时刻目标i在模型ri,k下的预测高斯伽玛混合分布为
Figure FDA0002400948370000022
预测存在概率为ρi,k|k-1(ri,k);其中,mi,k|k-1(ri,k)=Fk-1(ri,k)mi,k-1(ri,k-1),
Figure FDA0002400948370000023
Figure FDA0002400948370000024
γi,k|k-1(ri,k)=ργγi,k-1(ri,k-1),ηi,k|k-1(ri,k)=ρηηi,k-1(ri,k-1),ρi,k|k-1(ri,k)=ps,k(tk-tk-1)tk|k-1(ri,k|ri,k-1i,k-1(ri,k-1),Fk-1(ri,k)表示k-1时刻目标i在模型ri,k下的状态转移矩阵,Qk-1(ri,k)表示k-1时刻目标i在模型ri,k下的过程噪声协方差矩阵,上标T表示矩阵或向量的转置,ρα,ρβ,ργ,ρη表示传递因子,取值范围为(0,1],ps,k(tk-tk-1)表示目标的幸存概率,且
Figure FDA0002400948370000025
tk|k-1(ri,k|ri,k-1)表示模型间的马尔可夫转移概率,T表示采样周期,δ为预定的常数,l=1,...,ξ,ξ表示状态向量的维数;
步骤2、根据当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率,利用变分贝叶斯方法对当前时刻的测量数据进行序贯处理,得到当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率;
步骤3、对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率分别进行融合,得到当前时刻已存在目标的运动模型、高斯伽玛混合分布和存在概率;
步骤4、利用当前时刻的测量数据生成新生目标的高斯伽玛混合分布,为所述新生目标指定存在概率,并选取匀速运动模型作为新生目标的运动模型,将所述新生目标的运动模型、高斯伽玛混合分布及存在概率分别与所述当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,得到当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率;
步骤5、从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标,并将裁减后余下的目标的运动模型、高斯伽玛混合分布和存在概率作为滤波器下一次递归的输入;
从所述裁减后余下的目标中提取存在概率大于第二阈值的目标,并将提取出的目标的高斯伽玛混合分布作为当前时刻的输出,将输出的高斯伽玛混合分布中,高斯分布的均值作为当前时刻目标的状态估计。
2.如权利要求1所述的方法,其特征在于,设当前时刻的测量集为yk=(y1,k,...,yM,k),其中,M为当前时刻的测量数据的总数,则所述步骤2包括:
以当前时刻各个目标在Mr个运动模型下的预测高斯伽玛混合分布和预测存在概率,作为当前时刻各个目标在相应运动模型下的初始高斯伽玛混合分布和初始存在概率,则目标i在模型ri,k下的初始高斯伽玛混合分布取为
Figure FDA0002400948370000031
Figure FDA0002400948370000032
初始存在概率取为
Figure FDA0002400948370000033
其中i=1,2,...,Nk-1,ri,k=1,2,…,Mr
Figure FDA0002400948370000034
Figure FDA0002400948370000035
利用变分贝叶斯方法对第1个至第M个当前时刻的测量数据进行序贯处理,具体包括以下:
设第j个测量数据处理前,目标i在模型ri,k下的高斯伽玛混合分布及存在概率分别为
Figure FDA0002400948370000041
Figure FDA0002400948370000042
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr,j=1,2,…,M;
根据变分贝叶斯方法,计算得到第j个测量更新时,目标i在模型ri,k下的高斯伽玛混合分布为
Figure FDA0002400948370000043
存在概率分别为
Figure FDA0002400948370000044
其中,
Figure FDA0002400948370000045
Figure FDA0002400948370000046
Figure FDA0002400948370000047
表示伽玛函数,tr表示矩阵的迹,
Figure FDA0002400948370000048
表示均值向量,
Figure FDA0002400948370000049
表示协方差矩阵,
Figure FDA00024009483700000410
表示滤波器增益,伽玛分布的形状参数为
Figure FDA00024009483700000411
Figure FDA00024009483700000412
伽玛分布的尺度参数为
Figure FDA0002400948370000051
Figure FDA0002400948370000052
精度矩阵为
Figure FDA0002400948370000053
其中,
Figure FDA0002400948370000054
Figure FDA0002400948370000055
Figure FDA0002400948370000056
Hk为模型ri,k的观测矩阵,PD,k为目标的检测概率,λc,k为杂波密度,yj,k为当前时刻接收到的第j个测量数据,上标T表示矩阵或向量的转置,ξ为状态向量的维数;
Figure FDA0002400948370000057
则第j个测量数据处理后,目标i在模型ri,k下的的高斯伽玛混合分布为
Figure FDA0002400948370000058
Figure FDA0002400948370000059
目标i的存在概率为
Figure FDA00024009483700000510
其中
Figure FDA00024009483700000511
Figure FDA00024009483700000512
Figure FDA00024009483700000513
则第j个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布为
Figure FDA00024009483700000514
Figure FDA00024009483700000515
目标i的存在概率为
Figure FDA00024009483700000516
其中,
Figure FDA00024009483700000517
Figure FDA0002400948370000061
在第M个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布和存在概率分别为
Figure FDA0002400948370000062
Figure FDA0002400948370000063
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr
将第M个测量数据处理后,各个目标在不同模型下的高斯伽玛混合分布及存在概率分别作为当前时刻各个目标的更新高斯伽玛混合分布和存在概率,则所述当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布为
Figure FDA0002400948370000064
Figure FDA0002400948370000065
更新存在概率为
Figure FDA0002400948370000066
其中i=1,…,Nk-1,ri,k=1,2,…,Mr
Figure FDA0002400948370000067
Figure FDA0002400948370000068
Figure FDA0002400948370000069
3.根据权利要求2所述的方法,其特征在于,所述步骤3包括:
对当前时刻目标i在Mr个运动模型下的更新的高斯伽玛混合分布和存在概率分别进行融合,得到当前时刻已存在目标q的运动模型为rq,k、高斯伽玛混合分布为
Figure FDA00024009483700000610
和存在概率为ρq,k(rq,k),其中,q为当前时刻已存在目标的索引号,且q=1,…,Nk-1
Figure FDA00024009483700000611
Figure FDA00024009483700000612
伽玛分布的形状参数为
Figure FDA00024009483700000613
Figure FDA0002400948370000071
伽玛分布的尺度参数为
Figure FDA0002400948370000072
Figure FDA0002400948370000073
i=1,…,Nk-1,q=1,…,Nk-1
4.如权利要求3所述的方法,其特征在于,所述步骤4包括:
选取匀速运动模型作为新生目标的运动模型,设此时新生目标的运动模型为rj,k=1,其中,j=1,…,M;
利用M个测量数据生成当前时刻新生目标的高斯伽玛混合分布为
Figure FDA0002400948370000074
其中,j=1,…,M,
Figure FDA0002400948370000075
为第j个新生目标的协方差,
Figure FDA0002400948370000076
由第j个测量数据yj,k=[xj,k yj,k]T产生,且
Figure FDA0002400948370000077
Figure FDA0002400948370000078
Figure FDA0002400948370000079
为第j个新生目标分布函数中伽玛分布的形状参数,
Figure FDA00024009483700000710
Figure FDA00024009483700000711
为第j个新生目标分布函数中伽玛分布的尺度参数;
指定当前时刻各新生目标的存在概率为
Figure FDA00024009483700000712
其中,j=1,…,M,ρυ为所指定的存在概率;
将所述当前时刻已存在目标的运动模型与所述当前时刻新生目标的运动模型合并,得到当前时刻各个目标的运动模型为
Figure FDA00024009483700000713
其中,Nk=Nk-1+M;
将所述当前时刻已存在目标的高斯伽玛混合分布与所述当前时刻新生目标的高斯伽玛混合分布进行合并,得到当前时刻各个目标的高斯伽玛混合分布为
Figure FDA0002400948370000081
将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并,得到所述当前时刻各个目标的存在概率为
Figure FDA0002400948370000082
5.一种适用于闪烁噪声的多机动目标跟踪系统,其特征在于,所述系统包括:
预测模块,用于利用前一时刻各个目标的运动模型、以及高斯伽玛混合分布和存在概率产生伽玛分布的形状参数和尺度参数;根据前一时刻各个目标的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、模型间的马尔可夫转移概率、以及伽玛分布的形状参数和尺度参数,得到当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率;所述预测模块具体用于:
以k-1表示前一时刻,k表示当前时刻,tk-1表示前一时刻的时间,tk表示当前时刻的时间,ri,k-1表示前一时刻目标i的运动模型标识,且1≤ri,k-1≤Mr,ri,k表示当前时刻目标i的模型标识,且ri,k=1,2,…,Mr,其中,Mr表示运动模型的总数目;当前时刻的测量噪声服从ξ维的学生t分布,以S(yj,k;Hkxi,k,Ri,k,νi,k)表示当前时刻测量的概率密度函数,其中,yj,k表示当前时刻的一个测量,Hk表示测量矩阵,xi,k表示目标的状态,Ri,k表示精度矩阵,νi,k表示t分布的自由度,且
Figure FDA0002400948370000083
Figure FDA0002400948370000084
为精度矩阵Ri,k主对角线上的各个元素;
已知k-1时刻目标i的运动模型为ri,k-1,则在该模型ri,k-1下的高斯伽玛混合分布为
Figure FDA0002400948370000091
存在概率为ρi,k-1(ri,k-1);其中,N表示高斯分布,g表示伽玛分布,xi,k-1表示k-1时刻目标i的状态向量,mi,k-1(ri,k-1)和Pi,k-1(ri,k-1)分别表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中高斯分布的均值和协方差,
Figure FDA0002400948370000092
表示k-1时刻精度矩阵
Figure FDA0002400948370000093
的对角线元素,νi,k-1表示k-1时刻t分布的自由度,
Figure FDA0002400948370000094
和γi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的形状参数,
Figure FDA0002400948370000095
和ηi,k-1(ri,k-1)表示k-1时刻目标i在模型ri,k-1下的高斯伽玛混合分布中伽玛分布的尺度参数,ξ表示状态向量的维数,i=1,2,...,Nk-1,Nk-1表示前一时刻目标的总数;
根据k-1时刻目标i的运动模型、高斯伽玛混合分布和存在概率、当前时刻与前一时刻的时间差、以及模型间的马尔可夫转移概率,得到当前时刻目标i在模型ri,k下的预测高斯伽玛混合分布为
Figure FDA0002400948370000096
预测存在概率为ρi,k|k-1(ri,k);其中,mi,k|k-1(ri,k)=Fk-1(ri,k)mi,k-1(ri,k-1),
Figure FDA0002400948370000097
Figure FDA0002400948370000098
γi,k|k-1(ri,k)=ργγi,k-1(ri,k-1),ηi,k|k-1(ri,k)=ρηηi,k-1(ri,k-1),ρi,k|k-1(ri,k)=ps,k(tk-tk-1)tk|k-1(ri,k|ri,k-1i,k-1(ri,k-1),Fk-1(ri,k)表示k-1时刻目标i在模型ri,k下的状态转移矩阵,Qk-1(ri,k)表示k-1时刻目标i在模型ri,k下的过程噪声协方差矩阵,上标T表示矩阵或向量的转置,ρα,ρβ,ργ,ρη表示传递因子,取值范围为(0,1],ps,k(tk-tk-1)表示目标的幸存概率,且
Figure FDA0002400948370000099
tk|k-1(ri,k|ri,k-1)表示模型间的马尔可夫转移概率,T表示采样周期,δ为预定的常数,l=1,...,ξ,ξ表示状态向量的维数;
更新模块,用于根据当前时刻各个目标在不同运动模型下的预测高斯伽玛混合分布和预测存在概率,利用变分贝叶斯方法对当前时刻的测量数据进行序贯处理,得到当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率;
融合模块,用于对当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布和更新存在概率分别进行融合,得到当前时刻已存在目标的运动模型、高斯伽玛混合分布和存在概率;
生成模块,用于利用当前时刻的测量数据生成新生目标的高斯伽玛混合分布,为所述新生目标指定存在概率,并选取匀速运动模型作为新生目标的运动模型,将所述新生目标的运动模型、高斯伽玛混合分布及存在概率分别与所述当前时刻已存在目标的运动模型、高斯伽玛混合分布及存在概率进行合并,得到当前时刻各个目标的运动模型、高斯伽玛混合分布和存在概率;
提取模块,用于从所述当前时刻各个目标中裁减掉存在概率小于第一阈值的目标,并将裁减后余下的目标的运动模型、高斯伽玛混合分布和存在概率作为滤波器下一次递归的输入;从所述裁减后余下的目标中提取存在概率大于第二阈值的目标,并将提取出的目标的高斯伽玛混合分布作为当前时刻的输出,将输出的高斯伽玛混合分布中,高斯分布的均值作为当前时刻目标的状态估计。
6.如权利要求5所述的系统,其特征在于,设当前时刻的测量集为yk=(y1,k,...,yM,k),其中,M为当前时刻测量的总数,则所述更新模块具体用于:
以当前时刻各个目标在Mr个运动模型下的预测高斯伽玛混合分布和预测存在概率,作为当前时刻各个目标在相应运动模型下的初始高斯伽玛混合分布和初始存在概率,则目标i在模型ri,k下的初始高斯伽玛混合分布取为
Figure FDA0002400948370000101
Figure FDA0002400948370000102
初始存在概率取为
Figure FDA0002400948370000111
其中i=1,2,...,Nk-1,ri,k=1,2,…,Mr
Figure FDA0002400948370000112
Figure FDA0002400948370000113
利用变分贝叶斯方法对第1个至第M个当前时刻的测量数据进行序贯处理,具体包括以下:
设第j个测量数据处理前,目标i在模型ri,k下的高斯伽玛混合分布及存在概率分别为
Figure FDA0002400948370000114
Figure FDA0002400948370000115
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr,j=1,2,…,M;
根据变分贝叶斯方法,计算得到第j个测量更新时,目标i在模型ri,k下的高斯伽玛混合分布为
Figure FDA0002400948370000116
存在概率分别为
Figure FDA0002400948370000117
其中,
Figure FDA0002400948370000118
Figure FDA0002400948370000119
Figure FDA00024009483700001110
表示伽玛函数,tr表示矩阵的迹,
Figure FDA00024009483700001111
表示均值向量,
Figure FDA0002400948370000121
表示协方差矩阵,
Figure FDA0002400948370000122
表示滤波器增益,伽玛分布的形状参数为
Figure FDA0002400948370000123
Figure FDA0002400948370000124
伽玛分布的尺度参数为
Figure FDA0002400948370000125
Figure FDA0002400948370000126
精度矩阵为
Figure FDA0002400948370000127
其中,
Figure FDA0002400948370000128
Figure FDA0002400948370000129
Figure FDA00024009483700001210
Hk为模型ri,k的观测矩阵,PD,k为目标的检测概率,λc,k为杂波密度,yj,k为当前时刻接收到的第j个测量数据,上标T表示矩阵或向量的转置,ξ为状态向量的维数;
Figure FDA00024009483700001211
则第j个测量数据处理后,目标i在模型ri,k下的的高斯伽玛混合分布为
Figure FDA00024009483700001212
Figure FDA00024009483700001213
目标i的存在概率为
Figure FDA00024009483700001214
其中
Figure FDA00024009483700001215
Figure FDA00024009483700001216
Figure FDA00024009483700001217
则第j个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布为
Figure FDA0002400948370000131
Figure FDA0002400948370000132
目标i的存在概率为
Figure FDA0002400948370000133
其中,
Figure FDA0002400948370000134
Figure FDA0002400948370000135
在第M个测量数据处理后,目标i在模型ri,k下的高斯伽玛混合分布和存在概率分别为
Figure FDA0002400948370000136
Figure FDA0002400948370000137
其中,i=1,2,...,Nk-1,ri,k=1,2,…,Mr
将第M个测量数据处理后,各个目标在不同模型下的高斯伽玛混合分布及存在概率分别作为当前时刻各个目标的更新高斯伽玛混合分布和存在概率,则所述当前时刻各个目标在不同运动模型下的更新高斯伽玛混合分布为
Figure FDA0002400948370000138
Figure FDA0002400948370000139
更新存在概率为
Figure FDA00024009483700001310
其中i=1,…,Nk-1,ri,k=1,2,…,Mr
Figure FDA00024009483700001311
Figure FDA00024009483700001312
Figure FDA00024009483700001313
7.如权利要求6所述的系统,其特征在于,所述融合模块具体用于:
对当前时刻目标i在Mr个运动模型下的更新的高斯伽玛混合分布和存在概率分别进行融合,得到当前时刻已存在目标q的运动模型为rq,k、高斯伽玛混合分布为
Figure FDA00024009483700001314
和存在概率为ρq,k(rq,k),其中,q为当前时刻已存在目标的索引号,且q=1,…,Nk-1
Figure FDA0002400948370000141
Figure FDA0002400948370000142
伽玛分布的形状参数为
Figure FDA0002400948370000143
Figure FDA0002400948370000144
伽玛分布的尺度参数为
Figure FDA0002400948370000145
Figure FDA0002400948370000146
i=1,…,Nk-1,q=1,…,Nk-1
8.如权利要求7所述的系统,其特征在于,所述生成模块具体用于:
选取匀速运动模型作为新生目标的运动模型,设此时新生目标的运动模型为rj,k=1,其中,j=1,…,M;
利用M个测量数据生成当前时刻新生目标的高斯伽玛混合分布为
Figure FDA0002400948370000147
其中,j=1,…,M,
Figure FDA0002400948370000148
为第j个新生目标的协方差,
Figure FDA0002400948370000149
由第j个测量数据yj,k=[xj,k yj,k]T产生,且
Figure FDA00024009483700001410
Figure FDA00024009483700001411
Figure FDA00024009483700001412
为第j个新生目标分布函数中伽玛分布的形状参数,
Figure FDA00024009483700001413
Figure FDA00024009483700001414
为第j个新生目标分布函数中伽玛分布的尺度参数;
指定当前时刻各新生目标的存在概率为
Figure FDA00024009483700001415
其中,j=1,…,M,ρυ为所指定的存在概率;
将所述当前时刻已存在目标的运动模型与所述当前时刻新生目标的运动模型合并,得到当前时刻各个目标的运动模型为
Figure FDA00024009483700001416
其中,Nk=Nk-1+M;
将所述当前时刻已存在目标的高斯伽玛混合分布与所述当前时刻新生目标的高斯伽玛混合分布进行合并,得到当前时刻各个目标的高斯伽玛混合分布为
Figure FDA0002400948370000151
将所述当前时刻各个目标的存在概率与所述当前时刻新生目标的存在概率进行合并,得到所述当前时刻各个目标的存在概率为
Figure FDA0002400948370000152
CN201710804689.7A 2017-09-08 2017-09-08 一种适用于闪烁噪声的多机动目标跟踪方法及系统 Expired - Fee Related CN107462882B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201710804689.7A CN107462882B (zh) 2017-09-08 2017-09-08 一种适用于闪烁噪声的多机动目标跟踪方法及系统
PCT/CN2018/071580 WO2019047455A1 (zh) 2017-09-08 2018-01-05 一种适用于闪烁噪声的多机动目标跟踪方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710804689.7A CN107462882B (zh) 2017-09-08 2017-09-08 一种适用于闪烁噪声的多机动目标跟踪方法及系统

Publications (2)

Publication Number Publication Date
CN107462882A CN107462882A (zh) 2017-12-12
CN107462882B true CN107462882B (zh) 2020-06-02

Family

ID=60552223

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710804689.7A Expired - Fee Related CN107462882B (zh) 2017-09-08 2017-09-08 一种适用于闪烁噪声的多机动目标跟踪方法及系统

Country Status (2)

Country Link
CN (1) CN107462882B (zh)
WO (1) WO2019047455A1 (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107462882B (zh) * 2017-09-08 2020-06-02 深圳大学 一种适用于闪烁噪声的多机动目标跟踪方法及系统
CN109474892B (zh) * 2018-11-05 2020-07-10 浙江工商大学 基于信息形式的强鲁棒传感器网络目标跟踪方法
CN109657928B (zh) * 2018-11-28 2022-10-04 大连理工大学 一种车载传感器系统的闭环协同调度框架的协同调度方法
CN110390684A (zh) * 2019-07-16 2019-10-29 深圳大学 一种闪烁噪声下的多目标跟踪方法及系统
CN112733907A (zh) * 2020-12-31 2021-04-30 上海商汤临港智能科技有限公司 数据融合方法、装置、电子设备以及存储介质
CN113030940B (zh) * 2021-02-05 2023-11-17 河南科技大学 一种转弯机动下的多星凸型扩展目标跟踪方法
CN114445459B (zh) * 2022-01-16 2024-04-09 中国人民解放军空军工程大学 基于变分贝叶斯理论的连续-离散最大相关熵目标跟踪方法
CN114739383B (zh) * 2022-03-04 2024-05-14 中国人民解放军空军工程大学 一种基于多目标不确定性改进的gm-phd滤波方法
CN116500575B (zh) * 2023-05-11 2023-12-22 兰州理工大学 一种基于变分贝叶斯理论的扩展目标跟踪方法和装置

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9062978B2 (en) * 2011-05-31 2015-06-23 Massachusetts Institute Of Technology Tracking a body by nonlinear and non-Gaussian parametric filtering
CN103310115B (zh) * 2013-06-27 2015-12-23 西安电子科技大学 一种多目标跟踪的杂波估计方法
CN103345577B (zh) * 2013-06-27 2016-05-18 江南大学 变分贝叶斯概率假设密度多目标跟踪方法
CN103729637B (zh) * 2013-12-31 2017-01-11 西安工程大学 基于容积卡尔曼滤波的扩展目标概率假设密度滤波方法
CN104318059B (zh) * 2014-09-24 2018-02-02 深圳大学 用于非线性高斯系统的目标跟踪方法和跟踪系统
CN105719312B (zh) * 2016-01-19 2018-07-27 深圳大学 基于序贯贝叶斯滤波的多目标跟踪方法及跟踪系统
CN106772353B (zh) * 2016-11-29 2019-04-05 深圳大学 一种适用于闪烁噪声的多目标跟踪方法及系统
CN107462882B (zh) * 2017-09-08 2020-06-02 深圳大学 一种适用于闪烁噪声的多机动目标跟踪方法及系统

Also Published As

Publication number Publication date
WO2019047455A1 (zh) 2019-03-14
CN107462882A (zh) 2017-12-12

Similar Documents

Publication Publication Date Title
CN107462882B (zh) 一种适用于闪烁噪声的多机动目标跟踪方法及系统
CN109990786B (zh) 机动目标跟踪方法及装置
CN105719312B (zh) 基于序贯贝叶斯滤波的多目标跟踪方法及跟踪系统
WO2017124299A1 (zh) 基于序贯贝叶斯滤波的多目标跟踪方法及跟踪系统
Gao et al. Long short-term memory-based deep recurrent neural networks for target tracking
CN107677997B (zh) 基于GLMB滤波和Gibbs采样的扩展目标跟踪方法
CN104318059B (zh) 用于非线性高斯系统的目标跟踪方法和跟踪系统
CN106772353B (zh) 一种适用于闪烁噪声的多目标跟踪方法及系统
Aftab et al. A learning Gaussian process approach for maneuvering target tracking and smoothing
CN112613532B (zh) 基于雷达与循环神经网络补全红外融合的动目标跟踪方法
CN104156984A (zh) 一种不均匀杂波环境下多目标跟踪的概率假设密度方法
WO2018010099A1 (zh) 一种用于跟踪转弯机动目标的方法及其系统
CN103743401A (zh) 基于多模型航迹质量的异步融合方法
CN110111367A (zh) 模糊模型粒子滤波方法、装置、设备及存储介质
CN106168943A (zh) 一种用于跟踪转弯机动目标的方法及其系统
CN111830501B (zh) Hrrp历史特征辅助的信号模糊数据关联方法及系统
CN111798494A (zh) 广义相关熵准则下的机动目标鲁棒跟踪方法
Poddar et al. Tuning of GPS aided attitude estimation using evolutionary algorithms
Pollithy et al. Estimating uncertainties of recurrent neural networks in application to multitarget tracking
CN112415503A (zh) 一种基于目标重跟踪的多目标粒子滤波检测前跟踪方法
CN111811515A (zh) 一种基于高斯混合概率假设密度滤波器的多目标航迹提取方法
Liu et al. Multiple maneuvering target tracking by improved particle filter based on multiscan JPDA
Lan et al. Variational Nonlinear Kalman Filtering With Unknown Process Noise Covariance
CN116776158B (zh) 目标分类方法、装置及存储介质
CN116047495B (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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200602

Termination date: 20200908