CN115526323A - 低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法 - Google Patents

低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法 Download PDF

Info

Publication number
CN115526323A
CN115526323A CN202211106503.8A CN202211106503A CN115526323A CN 115526323 A CN115526323 A CN 115526323A CN 202211106503 A CN202211106503 A CN 202211106503A CN 115526323 A CN115526323 A CN 115526323A
Authority
CN
China
Prior art keywords
particle
maneuvering
target
particles
smoothing
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.)
Pending
Application number
CN202211106503.8A
Other languages
English (en)
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.)
Northwestern Polytechnical University
Shenzhen Institute of Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Shenzhen Institute of Northwestern Polytechnical 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 Northwestern Polytechnical University, Shenzhen Institute of Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202211106503.8A priority Critical patent/CN115526323A/zh
Publication of CN115526323A publication Critical patent/CN115526323A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • 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

Landscapes

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

Abstract

本发明涉及一种低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，针对多机动目标跟踪应用场景，建立具有未知混合状态的混杂系统模型，对多机动目标进行机动方式与运动状态的识别。设计出一种在粒子权重未退化时的低复杂度平滑计算方法，结合吉布斯采样、马尔可夫链蒙特卡洛方法估计未知的目标机动方式以及运动状态。与传统粒子平滑算法相比，本发明方法通过引入参考轨迹提高了目标机动方式和运动状态的估计精度，并通过低计算复杂度的粒子平滑算法提高了计算效率，因此该方法可在更短时间内达到相同精度的跟踪精度。

Description

低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法
技术领域
本发明属于混杂系统理论、随机过程分析在信号处理方面的应用，应用对象场景为目标机动模式识别与轨迹跟踪，涉及一种低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，具体方法涉及粒子平滑、吉布斯采样、马尔可夫链蒙特卡洛方法等。
背景技术
随着科技的迅猛发展，飞行器技术日益普及，给公共安全和国防安全带来巨大隐患和威胁。因此，迫切需要建立空中目标监视管理系统，而对上述飞行目标研究轨迹跟踪策略至关重要且势在必行。在这一领域中，由于目标运动方式具有不确定性和多变性，针对多机动目标的跟踪算法是实现目标精确、快速跟踪的基础，所以研究多机动目标轨迹跟踪算法具有重要的实际应用价值。
在相关的机动目标跟踪实际应用中，往往需要跟踪算法响应速度越快、跟踪精度越高，研究者们基于这两方面的需求提出了多种不同的算法。Lopez等(R.Lopez and P.Danès,"Low-Complexity IMM Smoothing for Jump Markov Nonlinear Systems,"in IEEETransactions on Aerospace and Electronic Systems,vol.53,no.3,pp.1261-1272,2017)提出了一种低复杂度的基于多交互模型的平滑估计算法，该算法能够以较快速度得到目标轨迹的估计结果。但是多交互模型算法对机动方式进行交互预测处理，使得该类算法无法精确判断目标机动方式，且对于高噪声、复杂机动等情形跟踪有精度较低、易发散等问题。另一种对于多机动目标跟踪的算法主要基于蒙特卡洛采样方法，通过建立具有跳变性质的非线性目标运动模型，利用粒子滤波或者粒子平滑算法对机动方式和目标轨迹进行估计。具有代表性的有传统的粒子平滑算法(T.T.Ashley and S. B.Andersson,"ASequential Monte Carlo framework for the system identification of jump Markovstate space models,"2014American Control Conference,2014,pp.1144-1149)、 Rao-Blackwellized粒子滤波与吉布斯先验采样结合的马尔可夫链蒙特卡洛方法(C. Cheng,J.-Y.Tourneret and X.Lu,"A Rao-Blackwellized Particle Filter With VariationalInference for State Estimation With Measurement Model Uncertainties,"in IEEEAccess, vol.8,pp.55665-55675,2020)等。
上述研究方法能够有效实现多机动目标跟踪，但是仍然存在如下问题：
(1)基于多目标交互模型的方法对机动方式转移概率较为敏感，导致机动方式估计不准确，对高噪声、机动情况复杂等情形无法保证跟踪精度，且易发散。
(2)基于蒙特卡洛采样方法虽然可以通过大量采样逼近真实的概率密度从而得到较为准确的跟踪结果，但是时间效率较低，常常达不到实时性要求。
发明内容
要解决的技术问题
为了避免现有技术的不足之处，本发明提出一种低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，针对基于传统粒子平滑算法的多机动目标跟踪算法粒子数少时跟踪精度低易发散、粒子数多时间效率差的问题，提出一种新的粒子平滑算法，实现相同跟踪精度下具有更高的时间效率。
技术方案
一种低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于步骤如下：
步骤1、基于混杂系统理论建立多机动目标运动模型：离散随机变量 R∈Ω＝{1,…,M}表示机动方式，其中M为机动方式种类的数量，连续随机变量X、Y 分别表示目标运动状态和观测量，其中R与X为未知，Y能够被观测到；n时刻的目标机动方式、运动状态、观测量的取值为rn、xn和yn，多机动目标运动模型为：
Figure BDA0003841864120000031
其中，alk为机动方式从l转移为k的转移概率，l,k∈Ω；非线性函数fn+1(xn,rn+1)为目标运动方程，由上一时刻运动状态和当前时刻的机动方式共同决定运动方程的具体形式；非线性函数gn+1(xn+1,rn+1)为观测方程；vx与vy分别是目标运动状态误差和观测误差，均为高斯白噪声；
系统模型参数为θ∈Θ＝{alk,θf|l,θg|l,Qx|l,Qy|l}l,k∈Ω，θf|l,θg|l分别为目标运动方程与观测方程的参数，Qx|l,Qy|l为误差vx、vy的协方差矩阵；连续随机变量X、Y分别表示目标运动状态和观测量
步骤2：设置粒子集合
Figure BDA0003841864120000032
其中K为粒子数量，粒子的数量代表了对目标运动状态和观测量连续随机变量X、Y的{x1:n,r1:n}概率分布的采样次数，数量越多表示采样得到的粒子集合越接近真实分布，
Figure BDA0003841864120000033
权重系数表示对应的粒子数值
Figure BDA0003841864120000034
与
Figure BDA0003841864120000035
出现的概率大小，第n时刻粒子权重满足
Figure BDA0003841864120000036
粒子重采样阈值hs，马尔可夫链蒙特卡洛法最大迭代次数λ，令迭代序号标志d＝1；
步骤3：设置目标运动状态和机动方式初始状态，根据多机动目标运动模型仿真得到观测量序列y1:N；
步骤4：采用粒子滤波获得马尔可夫链蒙特卡洛法初始参考轨迹{r′1:n,x′1:n}：
Figure BDA0003841864120000037
Figure BDA0003841864120000038
其中，
Figure BDA0003841864120000039
为满足
Figure BDA00038418641200000310
条件的狄拉克函数；通过引入参考轨迹{r′1:n,x′1:n}可使得估计得到的目标运动状态和机动方式更加接近于真实值；
步骤5：令αn(rn,xn)＝p(rn+1:N,yn+1:N|r1:n,x1:n,y1:n)，按照下式迭代计算所有时刻的αn(r′n,x′n)
Figure BDA0003841864120000041
步骤6：基于吉布斯采样从n＝2到n＝N执行滤波算法，更新粒子权重
Figure BDA0003841864120000042
步骤7：令
Figure BDA0003841864120000043
从时刻n＝N-1到n＝1迭代执行低计算复杂度的平滑算法，得到的粒子平滑的权重
Figure BDA0003841864120000044
步骤8：以步骤7得到的粒子平滑的权重
Figure BDA0003841864120000045
计算并更新所有n＝1,…N时刻的参考轨迹：
Figure BDA0003841864120000046
Figure BDA0003841864120000047
上式得到的机动方式参考轨迹r′n为使得由粒子平滑权重加权求和
Figure BDA0003841864120000048
最大值对应的机动方式；x′n为由粒子平滑权重加权求和
Figure BDA0003841864120000049
得到的加权求和平均值；
并将得到的参考轨迹记录为马尔可夫链蒙特卡洛方法第d次迭代的参考轨迹{r′1:n[d],x′1:n[d]}；令
Figure BDA00038418641200000410
若d≤λ，返回步骤5，否则执行步骤9；
步骤9：输出估计的机动方式和目标运动状态
Figure BDA00038418641200000411
Figure BDA00038418641200000412
所述非线性函数fn+1(xn,rn+1)的目标运动方程由上一时刻运动状态和当前时刻的机动方式共同决定运动方程的具体形式。
所述步骤2的粒子重采样阈值hs取值大于1小于5。
所述步骤4采用粒子滤波获得马尔可夫链蒙特卡洛法初始参考轨迹{r′1:n,x′1:n}的过程为：
步骤4.1：初始化所有K个粒子：令r1 (i)在Ω中随机取值，
Figure BDA00038418641200000413
为全零向量，
Figure BDA00038418641200000414
步骤4.2：按照下列方式迭代计算每个时刻的粒子
Figure BDA0003841864120000051
Figure BDA0003841864120000052
Figure BDA0003841864120000053
Figure BDA0003841864120000054
其中，符号‘～’表示左边数值为服从右边概率分布的随机采样，符号‘∝’表示左边数值正比于右边数值；
Figure BDA0003841864120000055
为根据上一时刻n-1的权重
Figure BDA0003841864120000056
的随机采样序号；函数q(·)表示概率密度采样函数，所以
Figure BDA0003841864120000057
表示按照机动方式从上一时刻为
Figure BDA0003841864120000058
变为
Figure BDA0003841864120000059
的概率进行采样，
Figure BDA00038418641200000510
表示在机动方式
Figure BDA00038418641200000511
条件下从上一时刻的目标运动状态
Figure BDA00038418641200000512
变为
Figure BDA00038418641200000513
进行采样；
Figure BDA00038418641200000514
表示观测值yn在
Figure BDA00038418641200000515
条件下的概率；若每一时刻迭代满足下式则对粒子按照权重
Figure BDA00038418641200000516
进行重新采样K个粒子，并令所有粒子权重
Figure BDA00038418641200000517
Figure BDA00038418641200000518
步骤4.3：按照下式计算初始参考轨迹{r′1:n,x′1:n}
Figure BDA00038418641200000519
Figure BDA00038418641200000520
其中，
Figure BDA00038418641200000521
为满足
Figure BDA00038418641200000522
条件的狄拉克函数；通过引入参考轨迹{r′1：n,x′1:n}可使得估计得到的目标运动状态和机动方式更加接近于真实值。
所述步骤6的基于吉布斯采样从n＝2到n＝N执行滤波算法，更新粒子权重
Figure BDA00038418641200000523
步骤为：
步骤6.1：按照下式采样K-1个粒子，
Figure BDA0003841864120000061
Figure BDA0003841864120000062
Figure BDA0003841864120000063
步骤6.2：依据参考轨迹{r′1:n,x′1:n}获得第K个粒子，首先计算先验权重：
Figure BDA0003841864120000064
Figure BDA0003841864120000065
根据
Figure BDA0003841864120000066
采样第K个粒子序号，并令
Figure BDA0003841864120000067
步骤6.3：按照下式更新粒子权重
Figure BDA0003841864120000068
Figure BDA0003841864120000069
步骤6.4：若
Figure BDA00038418641200000610
按照权重
Figure BDA00038418641200000611
重新采样K个粒子，并令所有粒子权重
Figure BDA00038418641200000612
返回步骤6.1直至n＝N。
所述步骤7计算复杂度的平滑算法：
步骤7.1：计算
Figure BDA00038418641200000613
步骤7.2：若
Figure BDA00038418641200000614
表示粒子平滑权重
Figure BDA00038418641200000615
出现了退化现象，意味着K个粒子中绝大部分粒子所代表的运动状态和机动方式数值出现的概率无限趋近于0，有效粒子数量过少；因此，为了保证计算准确性，按照传统粒子平滑算法计算粒子平滑的权重并归一化，如下：
Figure BDA00038418641200000616
若Neff＞hs，表示粒子平滑权重
Figure BDA00038418641200000617
没有退化，有效粒子足够多；因此，为了提高计算效率，使用低计算复杂度的粒子平滑方式计算粒子平滑的权重并归一化；
Figure BDA00038418641200000618
得到的粒子平滑的权重
Figure BDA00038418641200000619
可近似于看作对应粒子数值
Figure BDA00038418641200000620
和
Figure BDA00038418641200000621
在所有观测序列y1:N条件下出现的概率密度。
有益效果
本发明提出的一种低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，针对多机动目标跟踪应用场景，建立具有未知混合状态的混杂系统模型，对多机动目标进行机动方式与运动状态的识别。设计出一种在粒子权重未退化时的低复杂度平滑计算方法，结合吉布斯采样、马尔可夫链蒙特卡洛方法估计未知的目标机动方式以及运动状态。与传统粒子平滑算法相比，本发明方法通过引入参考轨迹提高了目标机动方式和运动状态的估计精度，并通过低计算复杂度的粒子平滑算法提高了计算效率，因此该方法可在更短时间内达到相同精度的跟踪精度。
附图说明
图1是方法流程图。
图2是机动方式估计误差随计算耗时的变化趋势比较图。
图中纵轴数据为通过设置不同粒子数量K和马尔可夫链蒙特卡洛方法迭代次数λ，并对各种不同K和λ取值执行跟踪算法100次，计算得到耗时在1秒内所有的机动方式估计误差的平均值。可以看到本发明方法相比传统算法在同样的计算耗时内可以更为准确地估计机动方式。
图3是运动状态估计误差随计算耗时的变化趋势比较图。
图中纵轴数据为通过设置不同粒子数量K和马尔可夫链蒙特卡洛方法迭代次数λ，并对各种不同K和λ取值执行跟踪算法100次，计算得到耗时在1秒内所有的运动状态估计误差的平均值。可以看到本发明方法相比传统算法在同样的计算耗时内可以更为准确地估计目标运动状态。
图4是计算耗时为49.5-50.5秒内轨迹跟踪结果比较图。
图中目标轨迹跟踪结果为当计算耗时为49.5-50.5秒的跟踪轨迹平均值。可以看到本发明方法能够得到的跟踪轨迹相比传统粒子平滑算法更贴近真实目标轨迹。
具体实施方式
现结合实施例、附图对本发明作进一步描述：
本发明提出基于低复杂度平滑估计的多机动目标跟踪方法，整个方法结合混杂系统理论、粒子平滑算法、马尔可夫链蒙特卡洛方法和吉布斯先验采样等，其具体步骤如下：
步骤一：
基于混杂系统理论建立多机动目标运动模型，令离散随机变量R∈Ω＝{1,…,M}表示机动方式，其中M为机动方式种类的数量，连续随机变量X、Y分别表示目标运动状态和观测量。在模型中，R与X未知，只有Y可被观测到，且在时刻n目标机动方式、运动状态、观测量的取值分别为rn、xn和yn，那么整个系统模型可以被表示为
Figure BDA0003841864120000081
其中，alk为机动方式从l转移为k的转移概率，l,k∈Ω；非线性函数fn+1(xn,rn+1)为目标运动方程，由上一时刻运动状态和当前时刻的机动方式共同决定运动方程的具体形式；非线性函数gn+1(xn+1,rn+1)为观测方程；
vx与vy分别是目标运动状态误差和观测误差，均为高斯白噪声。确立系统模型参数为θ∈Θ＝{alk,θf|l,θg|l,Qx|l,Qy|l}l,k∈Ω，θf|l,θg|l分别为目标运动方程与观测方程的参数，Qx|l,Qy|l为误差vx、vy的协方差矩阵。
步骤二：
设置粒子集合
Figure BDA0003841864120000082
其中K为粒子数量，粒子的数量代表了对{x1:n,r1:n} 概率分布的采样次数，数量越多表示采样得到的粒子集合越接近真实分布，
Figure BDA0003841864120000083
权重系数表示对应的粒子数值
Figure BDA0003841864120000084
与
Figure BDA0003841864120000085
出现的概率大小，第n时刻粒子权重满足
Figure BDA0003841864120000086
粒子重采样阈值hs(取值一般大于1小于5)，马尔可夫链蒙特卡洛法最大迭代次数λ，令迭代序号标志d＝1。
步骤三：
通过设置目标运动状态和机动方式初始状态，根据系统模型仿真得到观测量序列y1:N。
步骤四：
使用粒子滤波获得马尔可夫链蒙特卡洛法初始参考轨迹：
步骤4.1：初始化所有K个粒子：令r1 (i)在Ω中随机取值，
Figure BDA0003841864120000091
为全零向量，
Figure BDA0003841864120000092
步骤4.2：按照下列方式迭代计算每个时刻的粒子
Figure BDA0003841864120000093
Figure BDA0003841864120000094
Figure BDA0003841864120000095
Figure BDA0003841864120000096
其中，符号‘～’表示左边数值为服从右边概率分布的随机采样，符号‘∝’表示左边数值正比于右边数值；
Figure BDA0003841864120000097
为根据上一时刻n-1的权重
Figure BDA0003841864120000098
的随机采样序号；函数q(·)表示概率密度采样函数，所以
Figure BDA0003841864120000099
表示按照机动方式从上一时刻为
Figure BDA00038418641200000910
变为
Figure BDA00038418641200000911
的概率进行采样，
Figure BDA00038418641200000912
表示在机动方式
Figure BDA00038418641200000913
条件下从上一时刻的目标运动状态
Figure BDA00038418641200000914
变为
Figure BDA00038418641200000915
进行采样；
Figure BDA00038418641200000916
表示观测值yn在
Figure BDA00038418641200000917
条件下的概率。若每一时刻迭代满足下式则对粒子按照权重
Figure BDA00038418641200000918
进行重新采样K个粒子，并令所有粒子权重
Figure BDA00038418641200000919
Figure BDA00038418641200000920
步骤4.3：按照下式计算初始参考轨迹{r′1:n,x′1:n}
Figure BDA0003841864120000101
Figure BDA0003841864120000102
其中，
Figure BDA0003841864120000103
为满足
Figure BDA0003841864120000104
条件的狄拉克函数。通过引入参考轨迹{r′1:n,x′1:n}可使得估计得到的目标运动状态和机动方式更加接近于真实值。
步骤五：
令αn(rn,xn)＝p(rn+1:N,yn+1:N|r1:n,x1:n,y1:n)，按照下式迭代计算所有时刻的αn(r′n,x′n)
Figure BDA0003841864120000105
步骤六：
基于吉布斯采样从n＝2到n＝N执行滤波算法：
步骤6.1：按照下式采样K-1个粒子，
Figure BDA0003841864120000106
Figure BDA0003841864120000107
Figure BDA0003841864120000108
步骤6.2：依据参考轨迹{r′1:n,x′1:n}获得第K个粒子，首先计算先验权重：
Figure BDA0003841864120000109
Figure BDA00038418641200001010
根据
Figure BDA00038418641200001011
采样第K个粒子序号，并令
Figure BDA00038418641200001012
步骤6.3：按照下式更新粒子权重
Figure BDA00038418641200001013
Figure BDA00038418641200001014
步骤6.4：若
Figure BDA00038418641200001015
按照权重
Figure BDA00038418641200001016
重新采样K个粒子，并令所有粒子权重
Figure BDA00038418641200001017
返回步骤6.1直至n＝N。
步骤七：
令
Figure BDA00038418641200001018
从时刻n＝N-1到n＝1迭代执行低计算复杂度的平滑算法：
步骤7.1：计算
Figure BDA0003841864120000111
步骤7.2：若Neff＜hs，表示粒子平滑权重
Figure BDA0003841864120000112
出现了退化现象，意味着K个粒子中绝大部分粒子所代表的运动状态和机动方式数值出现的概率无限趋近于0，有效粒子数量过少。因此，为了保证计算准确性，按照传统粒子平滑算法计算粒子平滑的权重，如下：
Figure BDA0003841864120000113
并归一化。
若Neff＞hs，表示粒子平滑权重
Figure BDA0003841864120000114
没有退化，有效粒子足够多。因此，为了提高计算效率，使用低计算复杂度的粒子平滑方式计算粒子平滑的权重。
Figure BDA0003841864120000115
并归一化。
此步骤得到的粒子平滑的权重
Figure BDA0003841864120000116
可近似于看作对应粒子数值
Figure BDA0003841864120000117
和
Figure BDA0003841864120000118
在所有观测序列y1:N条件下出现的概率密度。
步骤八：
由步骤七得到的粒子平滑的权重
Figure BDA0003841864120000119
计算并更新所有n＝1,…N时刻的参考轨迹：
Figure BDA00038418641200001110
Figure BDA00038418641200001111
上式得到的机动方式参考轨迹r′n为使得由粒子平滑权重加权求和
Figure BDA00038418641200001112
最大值对应的机动方式；x′n为由粒子平滑权重加权求和
Figure BDA00038418641200001113
得到的加权求和平均值。并将得到的参考轨迹记录为马尔可夫链蒙特卡洛方法第d次迭代的参考轨迹{r′1:n[d],x′1:n[d]}。令
Figure BDA0003841864120000127
若d≤λ，返回步骤五，否则执行步骤九。
步骤九：
输出估计的机动方式和目标运动状态
Figure BDA0003841864120000121
Figure BDA0003841864120000122
现结合具体实例、附图对本发明作进一步描述：
令机动方式在有限集Ω∈{1,2,3}中选择，三种机动方式分别代表目标速度方向以10°/s角速率逆时针转动、直线运动(角速度为零)和以10°/s角速率顺时针转动，机动方式的转换具有马尔可夫性；目标运动状态表示为x＝[p1,υ1,p2,υ2]Τ，其中p1、p2分别表示目标位置在二维平面横轴和纵轴上的投影坐标，υ1和υ2分别表示目标速度在二维平面横轴和纵轴上的投影坐标；观测量二维向量，包含目标与坐标原点的距离，以及从坐标原点到目标位置的向量与坐标横轴间的夹角。那么可建立多机动目标在二维平面内的运动模型：
p(rn+1＝k|rn＝l)＝alk
Figure BDA0003841864120000123
Figure BDA0003841864120000124
Figure BDA0003841864120000125
Figure BDA0003841864120000126
其中ω为机动方式的角速度，T为观测时间间隔。
多机动目标跟踪的具体实现步骤如下：
步骤一：设置多机动目标跟踪系统参数如下表所示
Figure BDA0003841864120000131
如此可得粒子在给定上一时刻机动方式
Figure BDA0003841864120000132
条件下当前时刻机动方式
Figure BDA0003841864120000133
的采样
Figure BDA0003841864120000134
满足分布律
Figure BDA0003841864120000135
Figure BDA0003841864120000136
Figure BDA0003841864120000137
Figure BDA0003841864120000138
Figure BDA0003841864120000139
Figure BDA00038418641200001310
Figure BDA00038418641200001311
Figure BDA00038418641200001312
Figure BDA00038418641200001313
粒子在当前机动方式
Figure BDA00038418641200001314
和上一时刻运动状态
Figure BDA00038418641200001315
条件下的采样
Figure BDA00038418641200001316
满足以下高斯分布
Figure BDA00038418641200001317
粒子在当前机动方式
Figure BDA00038418641200001318
和运动状态
Figure BDA00038418641200001319
条件下的观测量概率密度为
Figure BDA0003841864120000141
步骤二：设粒子数量为K，粒子重采样阈值hs＝4，马尔可夫链蒙特卡洛法最大迭代次数λ，令迭代序号标志k＝1。
步骤三：令观测时间间隔T＝1s，目标初始运动状态x1＝[100,0,0,10]Τ，初始机动方式为r1＝1，仿真时间200秒共获得N＝200长度的数据{r1:N,x1:N,y1:N}。
步骤四：假设仅观测量y1:N已知，{r1:N,x1:N}未知，利用本发明的基于低计算复杂度粒子平滑的多机动目标跟踪算法估计目标的机动方式与运动状态，得到
Figure BDA0003841864120000142
并对比估计值
Figure BDA0003841864120000143
与真实值{r1:N,x1:N}之间的误差。

Claims (6)

1.一种低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于步骤如下：
步骤1、基于混杂系统理论建立多机动目标运动模型：离散随机变量R∈Ω＝{1,…,M}表示机动方式，其中M为机动方式种类的数量，连续随机变量X、Y分别表示目标运动状态和观测量，其中R与X为未知，Y能够被观测到；n时刻的目标机动方式、运动状态、观测量的取值为rn、xn和yn，多机动目标运动模型为：
Figure FDA0003841864110000011
其中，alk为机动方式从l转移为k的转移概率，l,k∈Ω；非线性函数fn+1(xn,rn+1)为目标运动方程，由上一时刻运动状态和当前时刻的机动方式共同决定运动方程的具体形式；非线性函数gn+1(xn+1,rn+1)为观测方程；vx与vy分别是目标运动状态误差和观测误差，均为高斯白噪声；
系统模型参数为θ∈Θ＝{alk,θf|l,θg|l,Qx|l,Qy|l}l,k∈Ω，θf|l,θg|l分别为目标运动方程与观测方程的参数，Qx|l,Qy|l为误差vx、vy的协方差矩阵；连续随机变量X、Y分别表示目标运动状态和观测量
步骤2：设置粒子集合
Figure FDA0003841864110000012
其中K为粒子数量，粒子的数量代表了对目标运动状态和观测量连续随机变量X、Y的{x1:n,r1:n}概率分布的采样次数，数量越多表示采样得到的粒子集合越接近真实分布，
Figure FDA0003841864110000013
权重系数表示对应的粒子数值
Figure FDA0003841864110000014
与
Figure FDA0003841864110000015
出现的概率大小，第n时刻粒子权重满足
Figure FDA0003841864110000016
粒子重采样阈值hs，马尔可夫链蒙特卡洛法最大迭代次数λ，令迭代序号标志d＝1；
步骤3：设置目标运动状态和机动方式初始状态，根据多机动目标运动模型仿真得到观测量序列y1:N；
步骤4：采用粒子滤波获得马尔可夫链蒙特卡洛法初始参考轨迹{r′1:n,x′1:n}：
Figure FDA0003841864110000021
Figure FDA0003841864110000022
其中，
Figure FDA0003841864110000023
为满足
Figure FDA0003841864110000024
条件的狄拉克函数；通过引入参考轨迹{r′1:n,x′1:n}可使得估计得到的目标运动状态和机动方式更加接近于真实值；
步骤5：令αn(rn,xn)＝p(rn+1:N,yn+1:N|r1:n,x1:n,y1:n)，按照下式迭代计算所有时刻的αn(r′n,x′n)
Figure FDA0003841864110000025
步骤6：基于吉布斯采样从n＝2到n＝N执行滤波算法，更新粒子权重
Figure FDA0003841864110000026
步骤7：令
Figure FDA0003841864110000027
从时刻n＝N-1到n＝1迭代执行低计算复杂度的平滑算法，得到的粒子平滑的权重
Figure FDA0003841864110000028
步骤8：以步骤7得到的粒子平滑的权重
Figure FDA0003841864110000029
计算并更新所有n＝1,…N时刻的参考轨迹：
Figure FDA00038418641100000210
Figure FDA00038418641100000211
上式得到的机动方式参考轨迹r′n为使得由粒子平滑权重加权求和
Figure FDA00038418641100000212
最大值对应的机动方式；x′n为由粒子平滑权重加权求和
Figure FDA00038418641100000213
得到的加权求和平均值；
并将得到的参考轨迹记录为马尔可夫链蒙特卡洛方法第d次迭代的参考轨迹{r′1:n[d],x′1:n[d]}；令
Figure FDA00038418641100000214
若d≤λ，返回步骤5，否则执行步骤9；
步骤9：输出估计的机动方式和目标运动状态
Figure FDA00038418641100000215
Figure FDA00038418641100000216
2.根据权利要求1所述低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于：所述非线性函数fn+1(xn,rn+1)的目标运动方程由上一时刻运动状态和当前时刻的机动方式共同决定运动方程的具体形式。
3.根据权利要求1所述低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于：所述步骤2的粒子重采样阈值hs取值大于1小于5。
4.根据权利要求1所述低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于：所述步骤4采用粒子滤波获得马尔可夫链蒙特卡洛法初始参考轨迹{r′1:n,x′1:n}的过程为：
步骤4.1：初始化所有K个粒子：令r1 (i)在Ω中随机取值，
Figure FDA0003841864110000031
为全零向量，
Figure FDA0003841864110000032
步骤4.2：按照下列方式迭代计算每个时刻的粒子
Figure FDA0003841864110000033
Figure FDA0003841864110000034
Figure FDA0003841864110000035
Figure FDA0003841864110000036
其中，符号‘～’表示左边数值为服从右边概率分布的随机采样，符号‘∝’表示左边数值正比于右边数值；
Figure FDA0003841864110000037
为根据上一时刻n-1的权重
Figure FDA0003841864110000038
的随机采样序号；函数q(·)表示概率密度采样函数，所以
Figure FDA0003841864110000039
表示按照机动方式从上一时刻为
Figure FDA00038418641100000310
变为
Figure FDA00038418641100000311
的概率进行采样，
Figure FDA00038418641100000312
表示在机动方式
Figure FDA00038418641100000313
条件下从上一时刻的目标运动状态
Figure FDA00038418641100000314
变为
Figure FDA00038418641100000315
进行采样；
Figure FDA00038418641100000316
表示观测值yn在
Figure FDA00038418641100000317
条件下的概率；若每一时刻迭代满足下式则对粒子按照权重
Figure FDA00038418641100000318
进行重新采样K个粒子，并令所有粒子权重
Figure FDA00038418641100000319
Figure FDA00038418641100000320
步骤4.3：按照下式计算初始参考轨迹{r′1:n,x′1:n}
Figure FDA0003841864110000041
Figure FDA0003841864110000042
其中，
Figure FDA0003841864110000043
为满足
Figure FDA0003841864110000044
条件的狄拉克函数；通过引入参考轨迹{r′1:n,x′1:n}可使得估计得到的目标运动状态和机动方式更加接近于真实值。
5.根据权利要求1所述低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于：所述步骤6的基于吉布斯采样从n＝2到n＝N执行滤波算法，更新粒子权重
Figure FDA0003841864110000045
步骤为：
步骤6.1：按照下式采样K-1个粒子，
Figure FDA0003841864110000046
Figure FDA0003841864110000047
Figure FDA0003841864110000048
步骤6.2：依据参考轨迹{r′1:n,x′1:n}获得第K个粒子，首先计算先验权重：
Figure FDA0003841864110000049
Figure FDA00038418641100000410
根据
Figure FDA00038418641100000411
采样第K个粒子序号，并令
Figure FDA00038418641100000412
步骤6.3：按照下式更新粒子权重
Figure FDA00038418641100000413
Figure FDA00038418641100000414
步骤6.4：若
Figure FDA00038418641100000415
按照权重
Figure FDA00038418641100000416
重新采样K个粒子，并令所有粒子权重
Figure FDA00038418641100000417
返回步骤6.1直至n＝N。
6.根据权利要求1所述低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法，其特征在于：所述步骤7计算复杂度的平滑算法：
步骤7.1：计算
Figure FDA00038418641100000418
步骤7.2：若Neff＜hs，表示粒子平滑权重
Figure FDA0003841864110000051
出现了退化现象，意味着K个粒子中绝大部分粒子所代表的运动状态和机动方式数值出现的概率无限趋近于0，有效粒子数量过少；按照传统粒子平滑算法计算粒子平滑的权重并归一化，如下：
Figure FDA0003841864110000052
若Neff＞hs，表示粒子平滑权重
Figure FDA0003841864110000053
没有退化，有效粒子足够多；使用低计算复杂度的粒子平滑方式计算粒子平滑的权重并归一化；
Figure FDA0003841864110000054
得到的粒子平滑的权重
Figure FDA0003841864110000055
近似于看作对应粒子数值
Figure FDA0003841864110000056
和
Figure FDA0003841864110000057
在所有观测序列y1:N条件下出现的概率密度。
CN202211106503.8A 2022-09-12 2022-09-12 低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法 Pending CN115526323A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211106503.8A CN115526323A (zh) 2022-09-12 2022-09-12 低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211106503.8A CN115526323A (zh) 2022-09-12 2022-09-12 低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法

Publications (1)

Publication Number Publication Date
CN115526323A true CN115526323A (zh) 2022-12-27

Family

ID=84698743

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211106503.8A Pending CN115526323A (zh) 2022-09-12 2022-09-12 低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法

Country Status (1)

Country Link
CN (1) CN115526323A (zh)

Similar Documents

Publication Publication Date Title
CN109508444B (zh) 区间量测下交互式多模广义标签多伯努利的快速跟踪方法
CN111178385B (zh) 一种鲁棒在线多传感器融合的目标跟踪方法
CN107402381B (zh) 一种迭代自适应的多机动目标跟踪方法
CN101975575B (zh) 基于粒子滤波的被动传感器多目标跟踪方法
CN105405151A (zh) 基于粒子滤波和加权Surf的抗遮挡目标跟踪方法
CN111711432B (zh) 一种基于ukf和pf混合滤波的目标跟踪算法
CN112819303B (zh) 基于pce代理模型的飞行器追踪效能评估方法及系统
CN106054167B (zh) 基于强度滤波器的多扩展目标跟踪方法
CN104484500A (zh) 一种基于拟合强化学习的空战行为建模方法
CN107797106A (zh) 一种加速em未知杂波估计的phd多目标跟踪平滑滤波方法
CN111562571A (zh) 一种未知新生强度的机动多目标跟踪与航迹维持方法
CN111798494A (zh) 广义相关熵准则下的机动目标鲁棒跟踪方法
CN111325776A (zh) 一种基于变分贝叶斯t分布卡尔曼滤波的phd多目标跟踪方法
Huang et al. A bank of maximum a posteriori estimators for single-sensor range-only target tracking
CN114236480A (zh) 一种机载平台传感器系统误差配准算法
CN115526323A (zh) 低计算复杂度粒子平滑估计的多机动目标轨迹跟踪方法
CN104880708B (zh) 一种可变数目机动目标跟踪方法
CN114415157B (zh) 一种基于水声传感器网络的水下目标多模型跟踪方法
CN111262556A (zh) 一种同时估计未知高斯测量噪声统计量的多目标跟踪方法
Zhang et al. Mobile robot localization based on gradient propagation particle filter network
CN115828533A (zh) 一种基于Student’s t分布的交互多模型鲁棒滤波方法
CN115578425A (zh) 一种应用在鱼苗计数器中基于无际卡尔曼滤波的动态追踪的方法
CN115544425A (zh) 一种基于目标信噪比特征估计的鲁棒多目标跟踪方法
CN104849697A (zh) 一种基于去偏坐标转换的α-β滤波方法
CN104467742A (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