CN115828533A - 一种基于Student’s t分布的交互多模型鲁棒滤波方法 - Google Patents

一种基于Student’s t分布的交互多模型鲁棒滤波方法 Download PDF

Info

Publication number
CN115828533A
CN115828533A CN202211422740.5A CN202211422740A CN115828533A CN 115828533 A CN115828533 A CN 115828533A CN 202211422740 A CN202211422740 A CN 202211422740A CN 115828533 A CN115828533 A CN 115828533A
Authority
CN
China
Prior art keywords
model
distribution
target
parameter
representing
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
CN202211422740.5A
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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202211422740.5A priority Critical patent/CN115828533A/zh
Publication of CN115828533A publication Critical patent/CN115828533A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • 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
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明属于智能信号处理技术领域，具体的说是一种基于Student's t分布的交互多模型鲁棒滤波方法。本发明是假设机动目标过程噪声与测量噪声均为拖尾噪声，通过将拖尾噪声建模为高斯伽马分布，并利用变分贝叶斯方法估计相关分布参数，最后对各个机动模型伪似然进行求解，从而解决了拖尾噪声背景下的机动目标多模型状态估计问题。该方法有适用范围广，鲁棒性强，估计精度高等特点，实现了拖尾噪声背景下的机动目标跟踪和参数估计，可以满足设计需求，具有良好的工程应用价值。

Description

一种基于Student’s t分布的交互多模型鲁棒滤波方法
技术领域
本发明属于智能信号处理技术领域，具体的说是一种基于Student's t分布的交互多模型鲁棒滤波方法。
背景技术
近年来，关于机动目标状态估计的方法已经得到的广泛的研究，针对目标状态估计方法，卡尔曼滤波由于其性能优越、易于实现和计算复杂度低，已成功地应用于许多实际应用，如目标跟踪、导航、定位等。对于高斯过程和测量噪声为线性状态空间模型时，卡尔曼滤波在最小均方误差(MMSE)方面是最优的。机动目标跟踪中的模型不确定性表现为目标的运动模型是时变的，不同模型之间的切换时间是未知的，由于模型的不确定性，卡尔曼滤波使用的模型可能与实际不一致，从而导致滤波逐渐发散。在模型不确定的情况下，多模型(MM)方法被认为是机动目标跟踪的主流方法，例如传统的交互多模型算法(IMM)、广义伪贝叶斯方法(GPB)等，这些方法建立有限的模型集来描述机动目标的不同运动模式，并通过各个模型的估计组合来实现最终的状态估计。
在一些工程应用中，如跟踪具有来自不可靠传感器的测量异常值的敏捷机动目标，重拖尾非高斯过程和测量噪声都存在，在重拖尾非高斯过程噪声和测量噪声的工程应用中，传统卡尔曼滤波器估计性能大幅恶化。为了解决重拖尾测量噪声线性系统的目标滤波问题，相关研究已经提出了大量的鲁棒滤波器，如Student's t混合滤波器、离群值鲁棒卡尔曼滤波器、Student's t和变分贝叶斯(VB)鲁棒卡尔曼滤波器，也有部分研究致力于解决量测拖尾噪声背景下的机动目标状态估计问题，但是目前关于过程噪声与量测噪声均为拖尾背景下机动目标状态估计问题仍亟待解决。
发明内容
针对上述问题，本发明提出一种基于Student's t分布的交互多模型鲁棒滤波方法。本发明假设机动目标状态估计中过程噪声与测量噪声均为拖尾噪声，通过将拖尾噪声Student's t分布近似建模为高斯-伽马分布，并利用变分贝叶斯方法估计相关分布参数，最后对各个机动模型伪似然进行求解，从而解决了拖尾噪声背景下的机动目标多模型状态估计问题。本发明方法解决实际工程应用中过程噪声与量测噪声均为拖尾噪声背景下的机动目标状态估计问题，具有较高的参数估计精度、状态估计精度、较好的对复杂应用背景的适应性和鲁棒性，可以满足工程应用需求。本发明大大提高了算法对复杂场景的适应性和鲁棒性，提高了拖尾噪声背景下的机动目标状态估计精度。
本发明的技术方案为：
一种基于Student's t分布的交互多模型鲁棒滤波方法，其特征在于，包括以下步骤：
S1、假设目标运动过程噪声与量测噪声均为拖尾噪声，初始化目标运动状态、量测状态的条件概率密度函数与其对应的相关参数(以下称为潜变量)先验分布，包括：
定义k时刻目标运动模型集合为rk∈{1,2,...,M}，其中M表示目标运动模型个数，定义k时刻运动状态为xk，则过程拖尾噪声背景下的第i个模型的运动状态条件概率密度函数表示为：
Figure BDA0003943165730000021
其中，z1:k-1表示目标从初始时刻到k-1时刻的量测集合向量，
Figure BDA0003943165730000022
表示k-1时刻模型i的状态估计值(由上一时刻获取)，fi(·)表示运动状态转移函数；
Figure BDA0003943165730000023
表示目标运动状态向量xk服从学生t(Student’s t)分布，其中
Figure BDA0003943165730000024
表示均值，
Figure BDA0003943165730000025
表示尺度矩阵，
Figure BDA0003943165730000026
表示自由度；
Figure BDA0003943165730000027
表示xk服从高斯分布，其中
Figure BDA0003943165730000028
表示高斯分布均值，
Figure BDA0003943165730000029
表示高斯分布协方差矩阵，其中
Figure BDA00039431657300000210
服从尺度矩阵
Figure BDA00039431657300000211
自由度为
Figure BDA00039431657300000212
的逆-wishart分布，即
Figure BDA00039431657300000213
G(x；a,b)表示伽马分布，其中a表示伽马分布形状参数，b表示伽马分布尺度参数，其中
Figure BDA00039431657300000214
服从形状参数为
Figure BDA00039431657300000215
尺度参数为
Figure BDA00039431657300000216
的伽马分布；∫·表示求积分操作；
量测噪声拖尾背景下的目标量测条件概率密度函数表示为：
Figure BDA00039431657300000217
其中，zk表示第k时刻的量测预测值，hi(·)表示第i个模型的量测转移函数，
Figure BDA0003943165730000031
表示zk服从Student’s-t分布，其中hi(xk)表示均值，
Figure BDA0003943165730000032
表示尺度矩阵，
Figure BDA0003943165730000033
表示自由度；
Figure BDA0003943165730000034
表示zk服从高斯分布，其中hi(xk)表示高斯分布均值，
Figure BDA0003943165730000035
表示高斯分布协方差矩阵，其中
Figure BDA0003943165730000036
服从尺度矩阵
Figure BDA0003943165730000037
自由度为
Figure BDA0003943165730000038
的逆-wishart分布，即
Figure BDA0003943165730000039
Figure BDA00039431657300000310
服从形状参数为
Figure BDA00039431657300000311
尺度参数为
Figure BDA00039431657300000312
的伽马分布。
S2、对目标各个模型的运动状态进行模型交互，包括：
结合下式计算目标模型j的预测概率
Figure BDA00039431657300000313
与模型交互概率
Figure BDA00039431657300000314
Figure BDA00039431657300000315
其中
Figure BDA00039431657300000316
表示k-1时刻模型i的概率更新值(由上一时刻获取)，
Figure BDA00039431657300000317
表示k时刻马尔科夫转移概率矩阵(TPM)；
结合下式计算目标第i个模型的运动状态交互结果
Figure BDA00039431657300000318
与状态协方差交互结果
Figure BDA00039431657300000319
Figure BDA00039431657300000320
其中，
Figure BDA00039431657300000321
与
Figure BDA00039431657300000322
分别表示k-1时刻模型j的状态估计值与状态估计协方差矩阵(由上一时刻获取)，∑·表示求和操作。
S3、对目标各个模型的运动状态值、运动状态协方差矩阵进行一步预测计算，包括：
结合下式计算目标第i个模型的运动状态一步预测值
Figure BDA00039431657300000323
与一步预测协方差矩阵
Figure BDA00039431657300000324
Figure BDA00039431657300000325
其中Fi表示模型运动状态转移函数fi(·)的雅克比矩阵，
Figure BDA00039431657300000326
表示k时刻模型i的过程噪声协方差矩阵，(·)T表示求转置操作。
S4、对目标各个模型潜变量进行模型交互，包括：
结合下式计算目标第i个模型的运动状态与量测潜变量逆-wishart分布Σk-1与Rk-1参数自由度交互结果
Figure BDA0003943165730000041
与尺度矩阵交互结果
Figure BDA0003943165730000042
Figure BDA0003943165730000043
其中
Figure BDA0003943165730000044
与
Figure BDA0003943165730000045
分别表示k-1时刻逆-wishart分布Σk-1与Rk-1的自由度，
Figure BDA0003943165730000046
与
Figure BDA0003943165730000047
分别表示k-1时刻逆-wishart分布Σk-1与Rk-1的尺度矩阵；
结合下式计算第i个模型的伽马分布
Figure BDA0003943165730000048
与
Figure BDA0003943165730000049
形状参数交互结果参数交互结果
Figure BDA00039431657300000410
与尺度参数交互结果
Figure BDA00039431657300000411
Figure BDA00039431657300000412
Figure BDA00039431657300000413
S5、对目标各个模型潜变量参数进行一步预测计算，包括：
结合下式计算第i个模型的逆-wishart分布
Figure BDA0003943165730000051
与
Figure BDA0003943165730000052
参数自由度预测值
Figure BDA0003943165730000053
与尺度矩阵预测值
Figure BDA0003943165730000054
Figure BDA0003943165730000055
其中nx表示运动状态向量维度，nz表示量测状态向量维度，ρ表示遗忘因子；
结合下式计算伽马分布
Figure BDA0003943165730000056
与
Figure BDA0003943165730000057
形状参数预测结果
Figure BDA0003943165730000058
与尺度参数预测结果
Figure BDA0003943165730000059
Figure BDA00039431657300000510
S6、结合变分贝叶斯算法求解目标状态更新方法与潜变量的超参数更新方法，
设定k时刻目标运动状态及潜变量集合Θk＝{xk,Σk,ξk,σk,Rk,λk,vk}，结合下式求解上述分布参数集合对应的联合概率密度函数，
Figure BDA00039431657300000511
其中z1:k表示目标从初始时刻到k时刻的量测值集合；
结合上述联合概率密度函数计算潜变量边缘概率密度(VB-marginal)并确定对应的超参数更新方法，
结合下式计算目标运动状态
Figure BDA00039431657300000512
边缘对数似然，
Figure BDA0003943165730000061
其中
Figure BDA0003943165730000062
E(·)表示求期望操作，Eq(x)(x)＝∫xq(x)dx，其中
Figure BDA0003943165730000063
Figure BDA0003943165730000064
表示求逆操作，
Figure BDA0003943165730000065
表示与
Figure BDA0003943165730000066
无关的常数项，结合上式可得目标第i个模型对应的状态更新值与状态更新协方差，
Figure BDA0003943165730000067
Figure BDA0003943165730000068
Figure BDA0003943165730000069
其中，Hi表示第i个模型量测转移函数的hi(·)的雅克比矩阵，
Figure BDA00039431657300000610
表示k时刻的新息协方差矩阵；
结合下式计算
Figure BDA00039431657300000611
边缘似然，
Figure BDA00039431657300000612
其中，det(·)表示求行列式操作，
Figure BDA00039431657300000613
表示联合概率密度函数中与
Figure BDA00039431657300000614
无关的期望常数项；其中，
Figure BDA00039431657300000615
综上可得到
Figure BDA00039431657300000616
对应的超参数更新步骤，
Figure BDA00039431657300000617
结合下式计算
Figure BDA00039431657300000618
边缘似然，
Figure BDA0003943165730000071
其中
Figure BDA0003943165730000072
Figure BDA0003943165730000073
表示与
Figure BDA0003943165730000074
无关的常数项；假定
Figure BDA0003943165730000075
确定ξk超参数更新步骤，
Figure BDA0003943165730000076
结合下式计算
Figure BDA0003943165730000077
边缘似然，
Figure BDA0003943165730000078
确定
Figure BDA0003943165730000079
超参数更新步骤，
Figure BDA00039431657300000710
其中
Figure BDA00039431657300000711
Ψ(·)表示双参数伽马分布。
获得量测相关潜变量的超参数更新步骤：
确定
Figure BDA00039431657300000712
超参数更新步骤，
Figure BDA00039431657300000713
其中，
Figure BDA00039431657300000714
假定
Figure BDA00039431657300000715
确定
Figure BDA00039431657300000716
超参数更新步骤，
Figure BDA0003943165730000081
其中
Figure BDA0003943165730000082
确定
Figure BDA0003943165730000083
超参数更新步骤:
Figure BDA0003943165730000084
其中
Figure BDA0003943165730000085
S7、计算目标各个运动模型对应的伪似然值log p(zk|z1:k-1)与k时刻对应的模型概率更新值
Figure BDA0003943165730000086
由变分贝叶斯原理结合下式获得伪似然结果
Figure BDA0003943165730000087
其中DKL(p||q)表示函数p与函数q之间的KL散度，结合S6对上式求解获得以下伪似然值，
Figure BDA0003943165730000088
根据获得的伪似然结果结合下式获得k时刻第i个模型的模型概率更新值
Figure BDA0003943165730000091
Figure BDA0003943165730000092
其中∝表示正比于。
S8、基于目标运动模型概率更新值xk|k与状态估计协方差矩阵Pk|k，
Figure BDA0003943165730000093
综上所述，得到目标在当前时刻的状态估计值与状态估计协方差矩阵，完成状态更新。
本发明的效益是，1.本发明在过程噪声与量测噪声均为拖尾分布情况下，基于变分贝叶斯技术对目标的过程噪声与量测噪声协方差分布参数进行实时估计，从而提高了拖尾噪声背景下的机动目标状态估计精度，并且提高了目标跟踪方法对复杂场景的适应性与鲁棒性；2.基于伪似然计算有效更新了机动目标运动模型概率，从而保证了机动目标跟踪的有效性。
附图说明
图1是本发明的整体流程图；
图2是采用本发明方法的目标状态真实值与量测值分布；
图3是采用最优状态估计、未进行参数估计与采用本发明方法时的目标距离估计精度，其中蒙特卡洛次数为500次；
图4是采用最优状态估计、未进行参数估计与采用本发明方法时的目标速度估计精度，其中蒙特卡洛次数为500次；
图5是采用最优状态估计、未进行参数估计与采用本发明方法时的模型概率估计值，其中蒙特卡洛次数为500次。
具体实施方式
下面结合附图和实施例，详细描述本发明的技术方案：
实施例
步骤1假设目标运动过程噪声与量测噪声均为拖尾噪声，初始化目标运动状态、量测状态的条件概率密度函数与其对应的相关参数(以下称为潜变量)先验分布，包括：
1.1.定义k时刻目标运动模型集合为rk∈{1,2,...,M}，其中M表示目标运动模型个数，定义k时刻运动状态为xk，则过程拖尾噪声背景下的第i个模型的运动状态条件概率密度函数表示为：
Figure BDA0003943165730000101
其中，z1:k-1表示目标从初始时刻到k-1时刻的量测集合向量，
Figure BDA0003943165730000102
表示k-1时刻模型i的状态估计值(由上一时刻获取)，fi(·)表示运动状态转移函数；
Figure BDA0003943165730000103
表示目标运动状态向量xk服从学生t(Student’s t)分布，其中
Figure BDA0003943165730000104
表示均值，
Figure BDA0003943165730000105
表示尺度矩阵，
Figure BDA0003943165730000106
表示自由度；
Figure BDA0003943165730000107
表示xk服从高斯分布，其中
Figure BDA0003943165730000108
表示高斯分布均值，
Figure BDA0003943165730000109
表示高斯分布协方差矩阵，其中
Figure BDA00039431657300001010
服从尺度矩阵
Figure BDA00039431657300001011
自由度为
Figure BDA00039431657300001012
的逆-wishart分布，即
Figure BDA00039431657300001013
G(x；a,b)表示伽马分布，其中a表示伽马分布形状参数，b表示伽马分布尺度参数，其中
Figure BDA00039431657300001014
服从形状参数为
Figure BDA00039431657300001015
尺度参数为
Figure BDA00039431657300001016
的伽马分布；∫·表示求积分操作；
本实例选用但不限于目标运动模型个数分别为匀速运动、已知转弯率的左右转弯模型，初始时刻目标运动各个模型的模型概率均初始化为以下形式：其中模型1表示目标为匀速运动，其状态转移矩阵与噪声扰动矩阵表示为：
Figure BDA00039431657300001017
T表示采样时间间隔，本实例选用但不限于T＝1s；模型2与3表示目标以已知转弯率的转弯模型，其状态转移矩阵表示为：
Figure BDA0003943165730000111
本实例选用但不限于w1＝0.1rad/s，w2＝-0.1rad/s；
本实例选用但不限于三种运动模型对应的状态向量分别为
Figure BDA0003943165730000112
其中x与y分别表示目标在x轴与y轴的距离，
Figure BDA00039431657300001119
与
Figure BDA00039431657300001120
分别表示目标在x轴与y轴的速度，状态转换对应的过程噪声方差均为1；
1.2量测噪声拖尾背景下的目标量测条件概率密度函数表示为：
Figure BDA0003943165730000113
其中，zk表示第k时刻的量测预测值，hi(·)表示第i个模型的量测转移函数，
Figure BDA0003943165730000114
表示zk服从学生t(Student’s-t)分布，其中hi(xk)表示均值，
Figure BDA0003943165730000115
表示尺度矩阵，
Figure BDA0003943165730000116
表示自由度；
Figure BDA0003943165730000117
表示zk服从高斯分布，其中hi(xk)表示高斯分布均值，
Figure BDA0003943165730000118
表示高斯分布协方差矩阵，其中
Figure BDA0003943165730000119
服从尺度矩阵
Figure BDA00039431657300001110
自由度为
Figure BDA00039431657300001111
的逆-wishart分布，即
Figure BDA00039431657300001112
Figure BDA00039431657300001113
服从形状参数为
Figure BDA00039431657300001114
尺度参数为
Figure BDA00039431657300001115
的伽马分布。
本实例选用但不限于量测转移函数为，
Figure BDA00039431657300001116
其中，sx与sy分别表示探测目标传感器的位置，量测噪声方差为[σR,σθ]＝[10m,0.2°]。
步骤2.对目标各个模型的运动状态及步骤1中的相关分布先验参数进行模型交互，包括：
2.1.结合下式计算目标模型j的预测概率
Figure BDA00039431657300001117
与模型交互概率
Figure BDA00039431657300001118
Figure BDA0003943165730000121
其中
Figure BDA0003943165730000122
表示k-1时刻模型i的概率(由上一时刻获取)，
Figure BDA0003943165730000123
表示k时刻马尔科夫转移概率矩阵(TPM)。本实例选用但不限于初始时刻各个模型初始概率
Figure BDA0003943165730000124
马尔科夫转移概率矩阵(TPM)具有以下形式：
Figure BDA0003943165730000125
2.2.结合步骤2.1获得的模型交互概率
Figure BDA0003943165730000126
与下式计算目标第i个模型的运动状态交互结果
Figure BDA0003943165730000127
与状态协方差交互结果
Figure BDA0003943165730000128
Figure BDA0003943165730000129
其中，
Figure BDA00039431657300001210
与
Figure BDA00039431657300001211
分别表示k-1时刻模型j的状态估计值与状态估计协方差矩阵(由上一时刻获取)，∑·表示求和操作。
步骤3.对目标各个模型的运动状态值、运动状态协方差矩阵进行一步预测计算，包括：
结合下式计算目标第i个模型的运动状态一步预测值
Figure BDA00039431657300001212
与一步预测协方差矩阵
Figure BDA00039431657300001213
Figure BDA00039431657300001214
其中Fi表示模型运动状态转移函数fi(·)的雅克比矩阵，
Figure BDA00039431657300001215
表示k时刻模型i的过程噪声协方差矩阵，(·)T表示求转置操作。
步骤4.对目标各个模型步骤1中的相关分布先验参数进行模型交互，包括：
4.1结合步骤2.1获得的模型交互概率
Figure BDA00039431657300001216
与下式计算目标第i个模型的运动状态与量测相关参数逆-wishart分布Σk-1与Rk-1参数自由度交互结果
Figure BDA00039431657300001217
与尺度矩阵交互结果
Figure BDA00039431657300001218
Figure BDA0003943165730000131
Figure BDA0003943165730000132
Figure BDA0003943165730000133
其中
Figure BDA0003943165730000134
与
Figure BDA0003943165730000135
分别表示k-1时刻逆-wishart分布Σk-1与Rk-1的自由度，
Figure BDA0003943165730000136
与
Figure BDA0003943165730000137
分别表示k-1时刻逆-wishart分布Σk-1与Rk-1的尺度矩阵；本实例对初始时刻尺度矩阵进行以下初始化，
Figure BDA0003943165730000138
其中nx表示运动状态向量维度，nz表示量测状态向量维度，ρ表示遗忘因子因子，τ表示调谐参数，本实用选用但不限于τ＝5,ρ＝0.98,nx＝4,nz＝2；
4.2.结合步骤2.1获得的模型交互概率
Figure BDA0003943165730000139
与下式计算第i个模型的伽马分布
Figure BDA00039431657300001310
与
Figure BDA00039431657300001311
形状参数交互结果参数交互结果
Figure BDA00039431657300001312
与尺度参数交互结果
Figure BDA00039431657300001313
Figure BDA00039431657300001314
Figure BDA00039431657300001315
Figure BDA0003943165730000141
本实例选用但不限于下述相关参数初始化结果
Figure BDA0003943165730000142
步骤5.对目标各个模型相关分布先验参数进行一步预测计算，包括：
5.1结合下式计算第i个模型的逆-wishart分布
Figure BDA0003943165730000143
与
Figure BDA00039431657300001413
参数自由度预测值
Figure BDA0003943165730000144
与尺度矩阵预测值
Figure BDA0003943165730000145
Figure BDA0003943165730000146
5.2结合下式计算伽马分布
Figure BDA0003943165730000147
与
Figure BDA0003943165730000148
形状参数预测结果
Figure BDA0003943165730000149
Figure BDA00039431657300001410
与尺度参数预测结果
Figure BDA00039431657300001411
Figure BDA00039431657300001412
步骤6.结合变分贝叶斯算法求解目标状态更新方法与相关参数分布的超参数更新方法，
6.1设定k时刻目标运动状态及分布参数集合Θk＝{xk,Σk,ξk,σk,Rk,λk,vk}，结合下式求解上述分布参数集合对应的联合概率密度函数，
Figure BDA0003943165730000151
其中z1:k表示目标从初始时刻到k时刻的量测值集合。
6.2结合上述联合概率密度函数计算目标分布参数边缘概率密度(VB-marginal)并确定对应的超参数更新方法，
6.2.1.结合下式计算目标运动状态
Figure BDA0003943165730000152
边缘对数似然，
Figure BDA0003943165730000153
其中
Figure BDA0003943165730000154
E(·)表示求期望操作，Eq(x)(x)＝∫xq(x)dx，其中
Figure BDA0003943165730000155
Figure BDA0003943165730000156
表示求逆操作，
Figure BDA0003943165730000157
表示与
Figure BDA00039431657300001512
无关的常数项，结合上式可得目标第i个模型对应的状态更新值与状态更新协方差，
Figure BDA0003943165730000158
Figure BDA0003943165730000159
Figure BDA00039431657300001510
其中，
Figure BDA00039431657300001511
表示k时刻的新息协方差矩阵，Hi表示第i个模型量测转移函数的hi(·)的雅克比矩阵，本实例选用以下形式，
Figure BDA0003943165730000161
6.2.2.结合下式计算
Figure BDA0003943165730000162
边缘似然，
Figure BDA0003943165730000163
其中，det(·)表示求行列式操作，
Figure BDA0003943165730000164
表示联合概率密度函数中与
Figure BDA0003943165730000165
无关的期望常数项，其中，
Figure BDA0003943165730000166
综上可得到
Figure BDA0003943165730000167
对应的超参数更新步骤，
Figure BDA0003943165730000168
Figure BDA0003943165730000169
6.2.3.结合下式计算
Figure BDA00039431657300001610
边缘似然，
Figure BDA00039431657300001611
其中
Figure BDA00039431657300001612
Figure BDA00039431657300001613
表示与
Figure BDA00039431657300001614
无关的常数项；假定
Figure BDA00039431657300001615
确定ξk超参数更新步骤，
Figure BDA00039431657300001616
6.2.4.结合下式计算
Figure BDA00039431657300001617
边缘似然，
Figure BDA0003943165730000171
确定
Figure BDA0003943165730000172
超参数更新步骤，
Figure BDA0003943165730000173
其中
Figure BDA0003943165730000174
Ψ(·)表示双参数伽马分布。
6.2.5.参考6.2.1～6.2.4可获得量测相关分布对应的超参数更新步骤：
确定
Figure BDA0003943165730000175
超参数更新步骤，
Figure BDA0003943165730000176
其中，
Figure BDA0003943165730000177
假定
Figure BDA0003943165730000178
确定
Figure BDA0003943165730000179
超参数更新步骤，
Figure BDA00039431657300001710
其中
Figure BDA00039431657300001711
确定
Figure BDA00039431657300001712
超参数更新步骤:
Figure BDA00039431657300001713
其中
Figure BDA00039431657300001714
步骤7.计算目标各个运动模型对应的伪似然值与k时刻对应的模型概率更新值，
7.1.结合步骤6结果与变分贝叶斯原理获得伪似然结果，
Figure BDA0003943165730000181
其中DKL(p||q)表示函数p与函数q之间的KL散度，结合步骤6对上式求解获得以下伪似然值，
Figure BDA0003943165730000182
7.2.根据步骤7.1获得的伪似然结果结合下式获得k时刻第i个模型的模型概率更新值
Figure BDA0003943165730000183
Figure BDA0003943165730000184
步骤8.结合下式目标运动模型概率更新值xk|k与状态估计协方差矩阵Pk|k，
Figure BDA0003943165730000185
Figure BDA0003943165730000186
综上所述，得到目标在当前时刻的状态估计值与状态估计协方差矩阵，完成状态估计。
通过下面的仿真示例证明本发明的实用性：
仿真条件及参数：
仿真场景为单目标跟踪场景，假设传感器位置为[sx,sy]＝[0,0]T，模型编号分别标记为1，2，3，模型转移过程为：{1,2,3,2,1,3,1}，每个模型持续时间为1→30s,2→40s,3→40s，过程噪声方差q＝1，量测噪声方差为[σR,σθ]＝[10m,0.2°]，总仿真时间长度为250s，过程噪声异常值比例为10％，对应的过程噪声方差为100；量测噪声异常值比例为10％，对应的量测噪声方差为40×[σR,σθ]。
仿真内容和结果分析：
图2是采用本发明方法中的目标运动轨迹变化图；
图3是不同方法下的距离均方根误差统计图；
均方根误差随时间变化曲线图，其中蒙特卡洛次数为500次，从图中可看出本发明方法中通过参数估计后相对于未进行参数估计方法距离估计误差均有明显的降低，证明了本发明方法的有效性。
图4是不同方法下的速度均方根误差统计图；
均方根误差随时间变化曲线图，其中蒙特卡洛次数为500次，从图中可看出本发明方法中通过参数估计后相对于未进行参数估计方法速度估计误差均有明显的降低，证明了本发明方法的有效性。
图5是采用本发明方法时模型概率估计均方根误差随时间变化曲线图，其中蒙特卡洛次数为500次，从图中可看出本发明方法具有较好的模型收敛概率。

Claims (1)

1.一种基于Student's t分布的交互多模型鲁棒滤波方法，其特征在于，包括以下步骤：
S1、定义目标运动过程噪声与量测噪声均为拖尾噪声，初始化目标运动状态、量测状态的条件概率密度函数与其对应的潜变量先验分布，包括：
定义k时刻目标运动模型集合为rk∈{1,2,...,M}，其中M表示目标运动模型个数，定义k时刻运动状态为xk，则过程拖尾噪声背景下的第i个模型的运动状态条件概率密度函数表示为：
Figure FDA0003943165720000011
其中，z1:k-1表示目标从初始时刻到k-1时刻的量测集合向量，
Figure FDA0003943165720000012
表示k-1时刻模型i的状态估计值，fi(·)表示运动状态转移函数；
Figure FDA0003943165720000013
表示目标运动状态向量xk服从学生Student’s-t分布，其中
Figure FDA0003943165720000014
表示均值，
Figure FDA0003943165720000015
表示尺度矩阵，
Figure FDA0003943165720000016
表示自由度；
Figure FDA0003943165720000017
表示xk服从高斯分布，其中
Figure FDA0003943165720000018
表示高斯分布均值，
Figure FDA0003943165720000019
表示高斯分布协方差矩阵，其中
Figure FDA00039431657200000110
服从尺度矩阵
Figure FDA00039431657200000111
自由度为
Figure FDA00039431657200000112
的逆-wishart分布，即
Figure FDA00039431657200000113
G(x；a,b)表示伽马分布，其中a表示伽马分布形状参数，b表示伽马分布尺度参数，其中
Figure FDA00039431657200000114
服从形状参数为
Figure FDA00039431657200000115
尺度参数为
Figure FDA00039431657200000116
的伽马分布；∫·表示求积分操作；
量测噪声拖尾背景下的目标量测条件概率密度函数为：
Figure FDA00039431657200000117
其中，zk表示第k时刻的量测预测值，hi(·)表示第i个模型的量测转移函数，
Figure FDA00039431657200000118
表示zk服从Student’s-t分布，其中hi(xk)表示均值，
Figure FDA00039431657200000119
表示尺度矩阵，
Figure FDA00039431657200000120
表示自由度；
Figure FDA00039431657200000121
表示zk服从高斯分布，其中hi(xk)表示高斯分布均值，
Figure FDA00039431657200000122
表示高斯分布协方差矩阵，其中
Figure FDA00039431657200000123
服从尺度矩阵
Figure FDA00039431657200000124
自由度为
Figure FDA00039431657200000125
的逆-wishart分布，即
Figure FDA00039431657200000126
Figure FDA00039431657200000127
服从形状参数为
Figure FDA00039431657200000128
尺度参数为
Figure FDA00039431657200000129
的伽马分布；
S2、对目标各个模型的运动状态进行模型交互，包括：
计算目标模型j的预测概率
Figure FDA0003943165720000021
与模型交互概率
Figure FDA0003943165720000022
Figure FDA0003943165720000023
其中
Figure FDA0003943165720000024
表示k-1时刻模型i的概率更新值，
Figure FDA0003943165720000025
表示k时刻马尔科夫转移概率矩阵；
计算目标第i个模型的运动状态交互结果
Figure FDA0003943165720000026
与状态协方差交互结果
Figure FDA0003943165720000027
Figure FDA0003943165720000028
其中，
Figure FDA0003943165720000029
与
Figure FDA00039431657200000210
分别表示k-1时刻模型j的状态估计值与状态估计协方差矩阵，∑·表示求和操作；
S3、对目标各个模型的运动状态值、运动状态协方差矩阵进行一步预测计算，包括：
结合下式计算目标第i个模型的运动状态一步预测值
Figure FDA00039431657200000211
与一步预测协方差矩阵
Figure FDA00039431657200000212
Figure FDA00039431657200000213
其中Fi表示模型运动状态转移函数fi(·)的雅克比矩阵，
Figure FDA00039431657200000214
表示k时刻模型i的过程噪声协方差矩阵，(·)T表示求转置操作；
S4、对目标各个模型潜变量进行模型交互，包括：
计算目标第i个模型的运动状态与量测潜变量逆-wishart分布Σk-1与Rk-1参数自由度交互结果
Figure FDA00039431657200000215
与尺度矩阵交互结果
Figure FDA00039431657200000216
Figure FDA0003943165720000031
其中
Figure FDA0003943165720000032
与
Figure FDA0003943165720000033
分别表示k-1时刻逆-wishart分布Σk-1与Rk-1的自由度，
Figure FDA0003943165720000034
与
Figure FDA0003943165720000035
分别表示k-1时刻逆-wishart分布Σk-1与Rk-1的尺度矩阵；
计算第i个模型的伽马分布
Figure FDA0003943165720000036
与
Figure FDA0003943165720000037
形状参数交互结果参数交互结果
Figure FDA0003943165720000038
与尺度参数交互结果
Figure FDA0003943165720000039
Figure FDA00039431657200000310
Figure FDA00039431657200000311
Figure FDA00039431657200000312
Figure FDA00039431657200000313
S5、对目标各个模型潜变量参数进行一步预测计算，包括：
计算第i个模型的逆-wishart分布
Figure FDA00039431657200000314
与
Figure FDA00039431657200000315
参数自由度预测值
Figure FDA00039431657200000316
与尺度矩阵预测值
Figure FDA00039431657200000317
Figure FDA0003943165720000041
其中nx表示运动状态向量维度，nz表示量测状态向量维度，ρ表示遗忘因子；
计算伽马分布
Figure FDA0003943165720000042
与
Figure FDA0003943165720000043
形状参数预测结果
Figure FDA0003943165720000044
与尺度参数预测结果
Figure FDA0003943165720000045
Figure FDA0003943165720000046
Figure FDA0003943165720000047
S6、结合变分贝叶斯算法求解目标状态更新方法与潜变量的超参数更新方法，
设定k时刻目标运动状态及潜变量集合Θk＝{xk,Σk,ξk,σk,Rk,λk,vk}，结合下式求解上述分布参数集合对应的联合概率密度函数，
Figure FDA0003943165720000048
其中z1:k表示目标从初始时刻到k时刻的量测值集合；
结合联合概率密度函数p(Θk|z1:k,rk＝i)计算潜变量边缘概率密度并确定对应的超参数更新方法；
计算目标运动状态
Figure FDA0003943165720000049
边缘对数似然，
Figure FDA00039431657200000410
其中
Figure FDA00039431657200000411
E(·)表示求期望操作，Eq(x)(x)＝∫xq(x)dx，其中
Figure FDA0003943165720000051
Figure FDA0003943165720000052
(·)-1表示求逆操作，
Figure FDA0003943165720000053
表示与
Figure FDA0003943165720000054
无关的常数项，结合上式可得目标第i个模型对应的状态更新值与状态更新协方差，
Figure FDA0003943165720000055
其中，Hi表示第i个模型量测转移函数的hi(·)的雅克比矩阵，
Figure FDA0003943165720000056
表示k时刻的新息协方差矩阵；
计算
Figure FDA0003943165720000057
边缘似然，
Figure FDA0003943165720000058
其中，det(·)表示求行列式操作，
Figure FDA0003943165720000059
表示联合概率密度函数中与
Figure FDA00039431657200000510
无关的期望常数项；其中，
Figure FDA00039431657200000511
综上可得到
Figure FDA00039431657200000512
对应的超参数更新步骤，
Figure FDA00039431657200000513
Figure FDA00039431657200000514
计算
Figure FDA00039431657200000515
边缘似然，
Figure FDA00039431657200000516
其中
Figure FDA00039431657200000517
Figure FDA00039431657200000518
表示与
Figure FDA00039431657200000519
无关的常数项；假定
Figure FDA00039431657200000520
确定ξk超参数更新步骤，
Figure FDA0003943165720000061
Figure FDA0003943165720000062
计算
Figure FDA0003943165720000063
边缘似然，
Figure FDA0003943165720000064
确定
Figure FDA0003943165720000065
超参数更新步骤，
Figure FDA0003943165720000066
其中
Figure FDA0003943165720000067
Ψ(·)表示双参数伽马分布；
获得量测相关潜变量的超参数更新步骤：
确定
Figure FDA0003943165720000068
超参数更新步骤，
Figure FDA0003943165720000069
其中，
Figure FDA00039431657200000610
假定
Figure FDA00039431657200000611
确定
Figure FDA00039431657200000612
超参数更新步骤，
Figure FDA00039431657200000613
其中
Figure FDA00039431657200000614
确定
Figure FDA00039431657200000615
超参数更新步骤:
Figure FDA0003943165720000071
其中
Figure FDA0003943165720000072
S7、计算目标各个运动模型对应的伪似然值logp(zk|z1:k-1)与k时刻对应的模型概率更新值
Figure FDA0003943165720000073
由变分贝叶斯原理结合下式获得伪似然结果
Figure FDA0003943165720000074
其中DKL(p||q)表示函数p与函数q之间的KL散度，结合S6对上式求解获得以下伪似然值，
Figure FDA0003943165720000075
根据获得的伪似然结果结合下式获得k时刻第i个模型的模型概率更新值
Figure FDA0003943165720000076
Figure FDA0003943165720000077
其中∝表示正比于；
S8、基于目标运动模型概率更新值xk|k与状态估计协方差矩阵Pk|k，
Figure FDA0003943165720000081
Figure FDA0003943165720000082
得到目标在当前时刻的状态估计值与状态估计协方差矩阵，完成状态更新。
CN202211422740.5A 2022-11-15 2022-11-15 一种基于Student’s t分布的交互多模型鲁棒滤波方法 Pending CN115828533A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211422740.5A CN115828533A (zh) 2022-11-15 2022-11-15 一种基于Student’s t分布的交互多模型鲁棒滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211422740.5A CN115828533A (zh) 2022-11-15 2022-11-15 一种基于Student’s t分布的交互多模型鲁棒滤波方法

Publications (1)

Publication Number Publication Date
CN115828533A true CN115828533A (zh) 2023-03-21

Family

ID=85528023

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211422740.5A Pending CN115828533A (zh) 2022-11-15 2022-11-15 一种基于Student’s t分布的交互多模型鲁棒滤波方法

Country Status (1)

Country Link
CN (1) CN115828533A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117109593A (zh) * 2023-10-20 2023-11-24 中国石油大学(华东) 一种基于鲁棒卡尔曼滤波的潜油机器人定位方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117109593A (zh) * 2023-10-20 2023-11-24 中国石油大学(华东) 一种基于鲁棒卡尔曼滤波的潜油机器人定位方法及系统
CN117109593B (zh) * 2023-10-20 2024-01-30 中国石油大学(华东) 一种基于鲁棒卡尔曼滤波的潜油机器人定位方法及系统

Similar Documents

Publication Publication Date Title
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN109990786B (zh) 机动目标跟踪方法及装置
CN111178385B (zh) 一种鲁棒在线多传感器融合的目标跟踪方法
CN110503071B (zh) 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法
CN107402381B (zh) 一种迭代自适应的多机动目标跟踪方法
CN103217175B (zh) 一种自适应容积卡尔曼滤波方法
CN106599368B (zh) 基于改进粒子提议分布和自适应粒子重采样的FastSLAM方法
CN111291471B (zh) 一种基于l1正则无迹变换的约束多模型滤波方法
Kim et al. Exactly Rao-Blackwellized unscented particle filters for SLAM
CN115828533A (zh) 一种基于Student’s t分布的交互多模型鲁棒滤波方法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
CN115204212A (zh) 一种基于stm-pmbm滤波算法的多目标跟踪方法
CN104777465B (zh) 基于b样条函数任意扩展目标形状及状态估计方法
CN111798494A (zh) 广义相关熵准则下的机动目标鲁棒跟踪方法
CN111291319A (zh) 一种应用于非高斯噪声环境下的移动机器人状态估计方法
CN112034445B (zh) 基于毫米波雷达的车辆运动轨迹跟踪方法和系统
CN111262556B (zh) 一种同时估计未知高斯测量噪声统计量的多目标跟踪方法
CN110912535B (zh) 一种新型无先导卡尔曼滤波方法
CN111340853B (zh) 基于ospa迭代的多传感器gmphd自适应融合方法
CN114061592B (zh) 基于多模型的自适应鲁棒auv导航方法
CN115905986A (zh) 一种基于联合策略的稳健卡尔曼滤波方法
CN111998854B (zh) 基于Cholesky分解计算的精确扩展Stirling插值滤波方法
CN113190960A (zh) 一种基于非等维状态混合估计的并行imm机动目标跟踪方法
CN115687890A (zh) 一种拖尾噪声背景下的机动目标分布式状态估计方法
Huang et al. Bearing-only target tracking using a bank of MAP estimators

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