CN114567288A - 基于变分贝叶斯的分布协同非线性系统状态估计方法 - Google Patents

基于变分贝叶斯的分布协同非线性系统状态估计方法 Download PDF

Info

Publication number
CN114567288A
CN114567288A CN202210088496.7A CN202210088496A CN114567288A CN 114567288 A CN114567288 A CN 114567288A CN 202210088496 A CN202210088496 A CN 202210088496A CN 114567288 A CN114567288 A CN 114567288A
Authority
CN
China
Prior art keywords
distribution
covariance
state
time
iteration
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
CN202210088496.7A
Other languages
English (en)
Other versions
CN114567288B (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.)
Henan University
Original Assignee
Henan 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 Henan University filed Critical Henan University
Priority to CN202210088496.7A priority Critical patent/CN114567288B/zh
Publication of CN114567288A publication Critical patent/CN114567288A/zh
Application granted granted Critical
Publication of CN114567288B publication Critical patent/CN114567288B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H03H17/0248—Filters characterised by a particular frequency response or filtering method
    • H03H17/0255—Filters based on statistics
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Physics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于变分贝叶斯的分布协同非线性系统状态估计方法，通过以变分贝叶斯为基础，针对分布协同非线性目标跟踪系统过程噪声时变和观测噪声具有随机异常值的情况，选择IW分布和student’t分布作为目标状态一步预测协方差的先验分布和量测的分布，通过定点迭代的方法在各局部滤波器中求解目标状态和噪声参数的近似后验分布，采用CI融合方法对各局部滤波器的目标状态估计加权融合得到全局最优估计，最后再将全局最优估计反馈给各局部滤波器，提高了滤波器的估计精度和稳定性。

Description

基于变分贝叶斯的分布协同非线性系统状态估计方法
技术领域
本发明涉及信号处理领域，具体涉及一种基于变分贝叶斯的分布协同非线性系统状态估计方法。
背景技术
无迹卡尔曼滤波算法(Unscented Kalman Filter,UKF)使用无迹变换来处理状态与协方差的非线性传递问题，用一系列确定sigma点来逼近状态的后验概率密度，具有很高的精度和稳定性，传统的UKF算法在系统噪声统计量已知且固定不变的非线性系统状态估计问题中能够取得较好的结果。但是在实际场景下，由于环境变化等原因，系统过程噪声统计量一般未知并且随着时间变化；由于外界扰动、传感器故障等原因，量测噪声异常值会随机出现；使用错误的噪声统计量和异常量测值对系统状态进行估计会导致估计结果变差，甚至引起滤波器发散。同时，受传感器精度影响，单传感器状态跟踪系统滤波精度有待提高。
在系统过程噪声时变、量测噪声存在随机异常值的情况下，要想对目标状态进行准确地估计，就要对时变的过程噪声和量测异常值进行估计，但是时变过程噪声与具有异常值的量测的概率分布难以直接得到，因此采用变分贝叶斯(Variational Bayesian,VB)近似方法来解决这个问题。通过对时变的过程噪声和异常量测值选取合适的先验分布，采用定点迭代的方法求取待估计的目标状态和噪声参数的近似后验分布，从而得到目标的状态估计。
针对单传感器状态跟踪系统滤波精度不足的问题，引入协方差交叉(Covarianceintersection,CI)融合算法，CI融合算法是一种经典的分布式融合算法，该算法能够对各局部传感器的估计结果进行加权融合得到全局最优估计，融合中心得到全局最优估计后，再将全局最优估计反馈给各局部滤波器，因此协方差交叉融合算法具有很高的估计精度，适合于传感器数量较少的分布式目标跟踪系统中。
发明内容
本发明的目的是提供一种基于变分贝叶斯的分布协同非线性系统状态估计方法，能够在系统过程噪声时变、量测噪声存在随机异常值情况下对目标状态进行准确地估计。
本发明采用的关键技术为：利用UKF处理非线性系统状态估计问题；利用变分贝叶斯方法解决了过程噪声时变和量测噪声具有随机异常值问题；利用定点迭代方法解决了目标状态和噪声参数相互耦合无法获取解析解的问题；利用CI融合方法解决了单传感器目标跟踪系统滤波精度不足问题。
一种基于变分贝叶斯的分布协同非线性系统状态估计方法，包括以下几个步骤：
1、建立分布协同非线性目标状态跟踪系统的动态模型；
所述步骤1具体包括以下步骤：
1.1、系统的动态空间模型为：
xk＝f(xk-1)+wk-1；
zk,l＝hl(xk)+vk,l l＝1,2,...s；
其中，k是离散时刻，xk、xk-1分别是k时刻和k-1时刻的状态向量，它们是n维的变量，zk,l是k时刻第l个传感器的量测向量，它是一个m维的变量，f(xk-1)是状态转移函数，hl(xk)是第l个传感器的量测函数；wk-1是k-1时刻到k时刻的零均值时变的过程噪声向量，它的期望协方差矩阵为Qk-1；vk,l是k时刻第l个传感器具有随机异常值的量测噪声向量，它的期望协方差矩阵为Rk,l；任意时刻的wk,vk,l和初值状态x0均互不相关。
2、在分布式非线性目标状态跟踪系统中，在k时刻，对于第l个局部滤波器(l＝1,2,…,s)，输入k-1时刻目标状态估计及其协方差及滤波器参数；
所述步骤2具体包括以下步骤：
2.1、输入：k-1时刻目标状态估计
Figure BDA0003488197370000031
及其对应的协方差Pk-1|k-1,l，k-1时刻的名义过程噪声协方差阵
Figure BDA0003488197370000032
k时刻该滤波器对应的量测噪声协方差阵Rk,l，k时刻该滤波器接收到的量测zk,l，调谐参数τl，student’t分布的自由度参数νl，以及变分迭代次数Nm；其中：k-1时刻的名义过程噪声协方差阵
Figure BDA0003488197370000033
是由于过程噪声统计量未知而根据经验选定k时刻的名义过程噪声协方差；τl的作用是协调模型先验信息和量测修正信息的权重。
3、在第l个局部滤波器中，通过UKF算法求k时刻目标状态一步预测
Figure BDA0003488197370000034
及其对应的协方差矩阵Pk|k-1,l；
所述步骤3具体包括以下步骤：
3-1、通过对第l个滤波器中k-1时刻目标状态估计
Figure BDA0003488197370000035
进行无迹变换产生2n+1个sigma点：
Figure BDA0003488197370000036
Figure BDA0003488197370000037
Figure BDA0003488197370000041
其中：n表示
Figure BDA0003488197370000042
的维度，
Figure BDA0003488197370000043
表示通过无迹变换产生的第j个sigma点，
Figure BDA0003488197370000044
表示第j个sigma点的权值，
Figure BDA0003488197370000045
表示第j个sigma点对应协方差阵的权值，
Figure BDA0003488197370000046
表示矩阵P方根的第j列，λ＝α2(n+κ)-n用来降低总的预测误差，α控制采样点的分布状态，κ的选取应保证(n+λ)P是半正定矩阵，一般取值0，β≥0合并方程中高阶项动差。
3-2、求状态预测及其对应的协方差阵：
Figure BDA0003488197370000047
Figure BDA0003488197370000048
Figure BDA0003488197370000049
其中：
Figure BDA00034881973700000410
是由于过程噪声统计量未知而根据经验选定的名义过程噪声协方差；
Figure BDA00034881973700000411
表示第j个sigma点的权值，
Figure BDA00034881973700000412
表示第j个sigma点对应协方差阵的权值。
4、选择逆威沙特(Inverse Wishart，IW)分布和student’t分布作为一步预测协方差的先验分布和量测的分布，并求取IW先验参数：
一步预测协方差的先验分布：
Figure BDA00034881973700000413
量测的分布：p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)
其中，
Figure BDA00034881973700000414
表示自由度参数为
Figure BDA00034881973700000415
和逆尺度参数为
Figure BDA00034881973700000416
的IW分布。St(zk,l；hl(xk),Rk,l,vl)表示均值为hl(xk)，尺度矩阵为Rk,l，自由度参数为v的student’t分布。
所述步骤4具体包括以下步骤：
4-1、选择IW分布作为一步预测协方差的先验分布：
Figure BDA0003488197370000051
其中，IW分布的自由度参数
Figure BDA0003488197370000052
逆尺度参数
Figure BDA0003488197370000053
通过如下方法得到：
对于一个IW分布：A～IW(A；t,T)，其期望可以写为：E[A-1]＝(t-n-1)Τ-1，t≥n+1，其中n是t的维度。因此，步骤3-2中的状态一步预测协方差Pk|k-1,l还可以表示为：
Figure BDA0003488197370000054
令：
Figure BDA0003488197370000055
则：
Figure BDA0003488197370000056
其中，nx是状态量x的维度，τl≥0是调谐参数，它的选取根据具体的情况而定。
4-2、使用student’t分布作为量测的分布：
p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)
St(zk,l；hl(xk),Rk,l,vl)表示均值为hl(xk)，尺度矩阵为Rk,l，自由度参数为v的student’t分布。各滤波器对目标状态估计独立进行，在各滤波器中，假设：
p(zk,l|xk,l)＝St(zk,l；hl(xk,l),Rk,l,vl)≈p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)，由于student’t分布的概率密度函数的闭合解难以求得，引入一个辅助随机变量λk,l，量测的概率密度函数可以写为：
Figure BDA0003488197370000057
其中，
Figure BDA0003488197370000061
表示形状参数为
Figure BDA0003488197370000062
和逆尺度参数为
Figure BDA0003488197370000063
的伽马分布。根据上一个式子，量测的概率密度函数最终可以表示为如下的分层高斯形式：
p(zk,l|xk,l,λk,l)＝N(zk,l；hl(xk,l),Rk,l/λk,l)
Figure BDA0003488197370000064
5、变分迭代初始化，根据步骤4得到IW先验参数tk|k-1,l,Tk|k-1,l；k时刻第l个局部滤波器中目标状态及协方差阵的迭代初值设置为：
Figure BDA0003488197370000065
student’t分布的辅助变量初值选取为
Figure BDA0003488197370000066
6、变分迭代：在k时刻、第i次变分迭代中(i＝1,2,…,Nm)，分别求取目标状态
Figure BDA0003488197370000067
一步预测协方差Pk|k-1,l和辅助变量λk,l的近似后验分布q(i)(xk|k,l),q(i)(Pk|k-1,l)和q(i)(λk,l)；
所述步骤6具体包括以下步骤：
6-1、固定k时刻第i-1次迭代状态估计的近似后验概率分布q(i-1)(xk,l)，将第i次迭代的状态一步预测协方差Pk|k-1,l的近似后验概率分布更新为
Figure BDA0003488197370000068
其中：自由度参数
Figure BDA0003488197370000069
逆尺度参数
Figure BDA00034881973700000610
k时刻第i次迭代一步预测协方差的逆矩阵期望为
Figure BDA00034881973700000611
第i次迭代一步预测协方差可以表示为：
Figure BDA00034881973700000612
6-2、固定k时刻第i-1次迭代状态估计的近似后验概率分布q(i-1)(xk,l)，将第i次迭代的辅助变量λk,l的近似后验概率分布更新为
Figure BDA00034881973700000613
其中：形状参数
Figure BDA0003488197370000071
逆尺度参数
Figure BDA0003488197370000072
Figure BDA0003488197370000073
Figure BDA0003488197370000074
是按照步骤3-1的方法对
Figure BDA0003488197370000075
进行无迹变换得到的第j个sigma点，
Figure BDA0003488197370000076
是对应的权值。
辅助变量期望
Figure BDA0003488197370000077
第i次迭代中辅助变量值表示为：
Figure BDA0003488197370000078
第i次迭代修正量测噪声协方差可以表示为：
Figure BDA0003488197370000079
6-3、固定k时刻第i次迭代状态一步预测协方差Pk|k-1,l的近似后验概率分布q(i)(Pk|k-1,l)、第i次迭代辅助变量λk,l的近似后验概率分布q(i)(λk,l)，将第i次迭代的对目标状态估计xk,l的近似后验概率分布更新为
Figure BDA00034881973700000710
其中：
Figure BDA00034881973700000711
分别表示第l个滤波器在k时刻第i次迭代中的目标状态估计及其对应的协方差。
按照步骤3-1的方法对状态一步预测
Figure BDA00034881973700000712
和第i次迭代得到的状态一步预测协方差
Figure BDA00034881973700000713
进行无迹变换重新获得状态一步预测sigma点以及量测预测
Figure BDA00034881973700000714
Figure BDA00034881973700000715
Figure BDA00034881973700000716
其中：n表示
Figure BDA00034881973700000717
的维度，
Figure BDA00034881973700000718
表示通过无迹变换产生的第j个sigma点，
Figure BDA00034881973700000719
表示矩阵P方根的第j列，
Figure BDA00034881973700000720
表示对第j个sigma点的量测预测，
Figure BDA00034881973700000721
表示第j个sigma点对应的权值，
Figure BDA00034881973700000722
表示第l个滤波器第k时刻的量测预测。
在UKF框架下求得
Figure BDA00034881973700000723
Figure BDA0003488197370000081
Figure BDA0003488197370000082
Figure BDA0003488197370000083
Figure BDA0003488197370000084
Figure BDA0003488197370000085
7、判断当前迭代次数i是否达到预设的最大变分迭代次数Nm，若是，则执行下一步；若不是，则i＝i+1，并重新执行步骤6；
8、对k时刻各局部滤波器的状态估计按照CI融合算法进行加权融合，再将融合结果反馈给各局部滤波器作为下一时刻的先验；
所述步骤8具体包括以下步骤：
8-1、局部滤波器状态估计的加权融合：
全局状态估计协方差：
Figure BDA0003488197370000086
全局状态估计：
Figure BDA0003488197370000087
8-2、将k时刻全局估计及其协方差按照一定的准则反馈给各局部滤波器：
状态估计反馈：
Figure BDA0003488197370000088
状态估计协方差反馈：
Figure BDA0003488197370000089
其中，αk,l为反馈权重系数，它们随着各个局部滤波器协方差的变化而变化，αk,l满足以下条件：
αk,1+αk,2+...+αk,l＝1
Figure BDA0003488197370000091
其中：||·||F表示表示Frobenius范数，即对于任意矩阵A：
Figure BDA0003488197370000092
9、输出第k时刻对目标状态的全局状态估计及其协方差：
Figure BDA0003488197370000093
10、判断是否达到预设仿真时长，若不是，则k＝k+1,重新执行步骤2；若是，结束执行；
本发明与现有技术相比，具有以下优点：
(1)引入变分贝叶斯方法，选择用IW分布和student’t分布分别作为目标状态一步预测协方差的先验分布和量测的分布，解决了过程噪声时变和量测噪声具有随机异常值问题；
(2)引入定点迭代方法，能够求得目标状态和噪声参数的近似后验概率分布，有效解决了目标状态和噪声参数相互耦合无法获取解析解的问题；
(3)引入CI融合算法，在分布协同目标跟踪系统中对所有局部滤波器的状态估计进行加权融合得到全局最优估计，有效解决了单传感器目标跟踪系统滤波精度不足问题；
附图说明
图1为本发明的流程图。
图2为本发明所述CI融合算法的结构图。
具体实施方式
下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
如图1流程图所示所示，本实施例提供一种基于变分贝叶斯的分布协同非线性系统状态估计方法，包括以下步骤：
1、建立分布协同非线性目标状态跟踪系统的动态模型；
2、在分布式非线性目标状态跟踪系统中，在k时刻，对于第l个局部滤波器(l＝1,2,…,s)，s为滤波器的个数，输入k-1时刻目标状态估计及其协方差及滤波器参数；
3、在第l个局部滤波器中通过UKF算法求目标状态一步预测
Figure BDA0003488197370000101
及其对应的协方差矩阵Pk|k-1,l；
4、选择逆威沙特(Inverse Wishart，IW)分布作为一步预测协方差的先验分布，并求取IW先验参数：
一步预测协方差的先验分布：
Figure BDA0003488197370000102
选择student’t分布作为量测的分布：
量测的分布：p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)
其中，
Figure BDA0003488197370000103
表示自由度参数为
Figure BDA0003488197370000104
和逆尺度参数为
Figure BDA0003488197370000105
的IW分布。St(zk,l；hl(xk),Rk,l,vl)表示均值为hl(xk)，尺度矩阵为Rk,l，自由度参数为v的student’t分布。
5、变分迭代初始化，根据步骤4得到IW先验参数tk|k-1,l,Tk|k-1,l；第l个局部滤波器中目标状态及协方差阵的迭代初值设置为：
Figure BDA0003488197370000106
student’t分布的辅助变量初值选取为
Figure BDA0003488197370000107
6、变分迭代：在k时刻、第i次变分迭代中(i＝1,2,…,Nm)，分别求取目标状态
Figure BDA0003488197370000108
的近似后验分布q(i)(xk|k,l)、一步预测协方差Pk|k-1,l的近似后验分布q(i)(Pk|k-1,l)和辅助变量λk,l的近似后验分布q(i)(λk,l)；
7、判断当前迭代次数i是否达到预设的次数Nm，若是，则执行下一步；若不是，则i＝i+1，并返回执行步骤6；
8、对各局部滤波器的状态估计按照CI融合算法进行加权融合，再将融合后的全局状态结果反馈给各局部滤波器作为下一时刻的先验；
9、输出第k时刻对目标状态的全局状态估计
Figure BDA0003488197370000111
及其协方差Pk|k,M；
10、判断是否达到预设仿真时长，若不是，则k＝k+1,重新执行步骤2；若是，结束执行；
其中，所述步骤1具体包括以下步骤：
1.1、系统的动态空间模型为：
xk＝f(xk-1)+wk-1；
zk,l＝hl(xk)+vk,l l＝1,2,...s；
其中，k是离散时刻，xk、xk-1分别是k时刻和k-1时刻的状态向量，它们是n维的变量，zk,l是k时刻第l个传感器的量测向量，它是一个m维的变量，f(xk-1)是状态转移函数，hl(xk)是第l个传感器的量测函数；wk-1是k-1时刻到k时刻的零均值时变的过程噪声向量，它的期望协方差矩阵为Qk-1；vk,l是k时刻第l个传感器具有随机异常值的量测噪声向量，它的期望协方差矩阵为Rk,l；任意时刻的wk,vk,l和初值状态x0均互不相关。
其中，所述步骤2具体包括以下步骤：
2.1、输入：k-1时刻目标状态估计
Figure BDA0003488197370000112
及其对应的协方差Pk-1|k-1,l，k-1时刻的名义过程噪声协方差阵
Figure BDA0003488197370000113
k时刻该滤波器对应的量测噪声协方差阵Rk,l，k时刻该滤波器接收到的量测zk,l，调谐参数τl，student’t分布的自由度参数νl，以及变分迭代次数Nm；其中：k-1时刻的名义过程噪声协方差阵
Figure BDA0003488197370000121
是由于过程噪声统计量未知而根据经验选定k时刻的名义过程噪声协方差；τl在算法中协调模型先验信息和量测修正信息的权重。
其中，所述步骤3具体包括以下步骤：
3-1、通过对第l个滤波器中k-1时刻目标状态估计
Figure BDA0003488197370000122
进行无迹变换产生2n+1个sigma点：
Figure BDA0003488197370000123
Figure BDA0003488197370000124
Figure BDA0003488197370000125
其中：n表示
Figure BDA0003488197370000126
的维度，
Figure BDA0003488197370000127
表示通过无迹变换产生的第j个sigma点，
Figure BDA0003488197370000128
表示第j个sigma点的权值，
Figure BDA0003488197370000129
表示第j个sigma点对应协方差阵的权值，
Figure BDA00034881973700001210
表示矩阵P方根的第j列，λ＝α2(n+κ)-n用来降低总的预测误差，α控制采样点的分布状态，κ的选取应保证(n+λ)P是半正定矩阵，一般取值0，β≥0合并方程中高阶项动差。
3-2、求状态预测及其对应的协方差阵：
Figure BDA00034881973700001211
Figure BDA00034881973700001212
Figure BDA0003488197370000131
其中：
Figure BDA0003488197370000132
是由于过程噪声统计量未知而根据经验选定的名义过程噪声协方差；
Figure BDA0003488197370000133
表示第j个sigma点的权值，
Figure BDA0003488197370000134
表示第j个sigma点对应协方差阵的权值。
其中，所述步骤4具体包括以下步骤：
4-1、选择IW分布作为一步预测协方差的先验分布：
Figure BDA0003488197370000135
其中，IW分布的自由度参数
Figure BDA0003488197370000136
逆尺度参数
Figure BDA0003488197370000137
通过如下方法得到：
对于对于一个IW分布：A～IW(A；t,T)，其期望可以写为：E[A-1]＝(t-n-1)Τ-1，t≥n+1，其中n是t的维度。因此，步骤3-2中的状态一步预测协方差Pk|k-1,l还可以表示为：
Figure BDA0003488197370000138
令：
Figure BDA0003488197370000139
则：
Figure BDA00034881973700001310
其中，nx是状态量x的维度，τl≥0是调谐参数，它的选取根据具体的情况而定。
4-2、使用student’t分布作为量测的分布：
p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)
St(zk,l；hl(xk),Rk,l,vl)表示均值为hl(xk)，尺度矩阵为Rk,l，自由度参数为v的student’t分布。各滤波器对目标状态估计独立进行，在各滤波器中，假设：
p(zk,l|xk,l)＝St(zk,l；hl(xk,l),Rk,l,vl)≈p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)，由于student’t分布的概率密度函数的闭合解难以求得，引入一个辅助随机变量λk,l，量测的概率密度函数可以写为：
Figure BDA0003488197370000141
其中，
Figure BDA0003488197370000142
表示形状参数为
Figure BDA0003488197370000143
和逆尺度参数为
Figure BDA0003488197370000144
的伽马分布。α和β分别为伽马分布的形状参数和逆尺度参数。根据上一个式子，量测的概率密度函数最终可以表示为如下的分层高斯形式：
p(zk,l|xk,l,λk,l)＝N(zk,l；hl(xk,l),Rk,l/λk,l)
Figure BDA0003488197370000145
其中，所述步骤6具体包括以下步骤：
6-1、固定k时刻第i次迭代状态估计的近似后验概率分布q(i-1)(xk,l)，将第i次迭代的状态一步预测协方差Pk|k-1,l的近似后验概率分布更新为
Figure BDA0003488197370000146
其中：自由度参数
Figure BDA0003488197370000147
逆尺度参数
Figure BDA0003488197370000148
k时刻第i次迭代一步预测协方差的逆矩阵期望为
Figure BDA0003488197370000149
第i次迭代一步预测协方差可以表示为：
Figure BDA00034881973700001410
6-2、固定k时刻第i-1次迭代状态估计的近似后验概率分布q(i-1)(xk,l)，将第i次迭代的辅助变量λk,l的近似后验概率分布更新为
Figure BDA00034881973700001411
其中：形状参数
Figure BDA00034881973700001412
逆尺度参数
Figure BDA00034881973700001413
Figure BDA00034881973700001414
Figure BDA00034881973700001415
是按照步骤3-1的方法对
Figure BDA00034881973700001416
进行无迹变换得到的第j个sigma点，
Figure BDA00034881973700001417
是对应的权值。
辅助变量期望
Figure BDA0003488197370000151
第i次迭代中辅助变量值表示为：
Figure BDA0003488197370000152
第i次迭代修正量测噪声协方差可以表示为：
Figure BDA0003488197370000153
6-3、固定k时刻第i次迭代状态一步预测协方差Pk|k-1,l的近似后验概率分布q(i)(Pk|k-1,l)、第i次迭代辅助变量λk,l的近似后验概率分布q(i)(λk,l)，将第i次迭代的对目标状态估计xk,l的近似后验概率分布更新为
Figure BDA0003488197370000154
其中：
Figure BDA0003488197370000155
分别表示第l个滤波器在k时刻第i次迭代中的目标状态估计及其对应的协方差。
按照步骤3-1的方法对状态一步预测
Figure BDA0003488197370000156
和第i次迭代得到的状态一步预测协方差
Figure BDA0003488197370000157
进行无迹变换重新获得状态一步预测sigma点以及量测预测
Figure BDA0003488197370000158
Figure BDA0003488197370000159
Figure BDA00034881973700001510
其中：n表示
Figure BDA00034881973700001511
的维度，
Figure BDA00034881973700001512
表示通过无迹变换产生的第j个sigma点，
Figure BDA00034881973700001513
表示矩阵P方根的第j列，
Figure BDA00034881973700001514
表示对第j个sigma点的量测预测，
Figure BDA00034881973700001515
表示第j个sigma点对应的权值，
Figure BDA00034881973700001516
表示第l个滤波器第k时刻的量测预测。
在UKF框架下求得
Figure BDA00034881973700001517
Figure BDA00034881973700001518
Figure BDA00034881973700001519
Figure BDA00034881973700001520
Figure BDA0003488197370000161
Figure BDA0003488197370000162
其中，所述步骤8具体包括以下步骤：
8-1、k时刻局部滤波器状态估计的加权融合：
全局状态估计协方差：
Figure BDA0003488197370000163
全局状态估计：
Figure BDA0003488197370000164
8-2、将k时刻全局估计及其协方差按照一定的准则反馈给各局部滤波器：
状态估计反馈：
Figure BDA0003488197370000165
状态估计协方差反馈：
Figure BDA0003488197370000166
其中，αk,l为反馈权重系数，它们随着各个局部滤波器协方差的变化而变化，αk,l满足以下条件：
αk,1+αk,2+...+αk,l＝1
Figure BDA0003488197370000167
其中：||·||F表示表示Frobenius范数，即对于任意矩阵A：
Figure BDA0003488197370000168
本发明以变分贝叶斯为基础，针对分布协同非线性目标跟踪系统过程噪声时变和观测噪声具有随机异常值的情况，选择IW分布和student’t分布作为目标状态一步预测协方差的先验分布和量测的分布，通过定点迭代的方法在各局部滤波器中求解目标状态和噪声参数的近似后验分布，采用CI融合方法对各局部滤波器的目标状态估计加权融合得到全局最优估计，最后再将全局最优估计反馈给各局部滤波器，提高了滤波器的估计精度和稳定性。

Claims (7)

1.一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：包括以下步骤：
1)、建立分布协同非线性目标状态跟踪系统的动态空间模型；
2)、在分布式非线性目标状态跟踪系统中，在k时刻，对于第l个局部滤波器，其中l＝1,2,…,s，s是滤波器的个数，输入k-1时刻目标状态估计及其协方差及滤波器参数；
3)、在第l个局部滤波器中通过UKF算法求目标状态一步预测
Figure FDA0003488197360000011
及其对应的协方差矩阵Pk|k-1,l；
4)、选择逆威沙特分布和student’t分布作为一步预测协方差的先验分布和量测的分布，并求取逆威沙特先验参数：
一步预测协方差的先验分布：
Figure FDA0003488197360000012
量测的分布：p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)
其中，
Figure FDA0003488197360000013
表示自由度参数为
Figure FDA0003488197360000014
和逆尺度参数为
Figure FDA0003488197360000015
的IW分布；St(zk,l；hl(xk),Rk,l,vl)表示均值为hl(xk)，尺度矩阵为Rk,l，自由度参数为v的student’t分布；
5)、变分迭代初始化，根据步骤4得到IW先验参数tk|k-1,l,Tk|k-1,l；第l个局部滤波器中目标状态及协方差阵的迭代初值设置为：
Figure FDA0003488197360000016
student’t分布的辅助变量初值选取为
Figure FDA0003488197360000017
6)、变分迭代：在k时刻、第i次变分迭代中(i＝1,2,…,Nm)，分别求取目标状态
Figure FDA0003488197360000018
一步预测协方差Pk|k-1,l和辅助变量λk,l的近似后验分布q(i)(xk|k,l),q(i)(Pk|k-1,l)和q(i)(λk,l)；
7)、判断当前迭代次数i是否达到Nm，若是，则执行下一步；若不是，则i＝i+1，并重新执行步骤6；
8)、对各局部滤波器的状态估计按照CI融合算法进行加权融合，再将融合结果反馈给各局部滤波器作为下一时刻的先验；
9)、输出第k时刻对目标状态的全局状态估计及其协方差：
Figure FDA0003488197360000019
10)、判断是否达到预设仿真时长，若不是，则k＝k+1,重新执行步骤2；若是，结束执行。
2.根据权利要求1所述的一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：所述步骤1)中所述系统的动态空间模型为：
xk＝f(xk-1)+wk-1；
zk,l＝hl(xk)+vk,l l＝1,2,...s；
其中，k是离散时刻，xk、xk-1分别是k时刻和k-1时刻的状态向量，它们是n维的变量，zk,l是k时刻第l个传感器的量测向量，它是一个m维的变量，f(xk-1)是状态转移函数，hl(xk)是第l个传感器的量测函数；wk-1是k-1时刻到k时刻的零均值时变的过程噪声向量，它的期望协方差矩阵为Qk-1；vk,l是k时刻第l个传感器具有随机异常值的量测噪声向量，它的期望协方差矩阵为Rk,l；任意时刻的wk,vk,l和初值状态x0均互不相关。
3.根据权利要求1所述的一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：所述步骤2)具体包括以下步骤：
输入：k-1时刻目标状态估计
Figure FDA0003488197360000021
及其对应的协方差Pk-1|k-1,l，k-1时刻的名义过程噪声协方差阵
Figure FDA0003488197360000022
k时刻该滤波器对应的量测噪声协方差阵Rk,l，k时刻该滤波器接收到的量测zk,l，调谐参数τl，student’t分布的自由度参数νl，以及变分迭代次数Nm；其中：k-1时刻的名义过程噪声协方差阵
Figure FDA0003488197360000023
是由于过程噪声统计量未知而根据经验选定k时刻的名义过程噪声协方差；τl的作用是协调模型先验信息和量测修正信息的权重。
4.根据权利要求1所述的一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：所述步骤3)具体包括以下步骤：
3-1)、对第l个滤波器中k-1时刻目标状态估计
Figure FDA0003488197360000024
通过无迹变换产生2n+1个sigma点：
Figure FDA0003488197360000025
Figure FDA0003488197360000031
Figure FDA0003488197360000032
其中：n表示
Figure FDA0003488197360000033
的维度，
Figure FDA0003488197360000034
表示通过无迹变换产生的第j个sigma点，
Figure FDA0003488197360000035
表示第j个sigma点的权值，
Figure FDA0003488197360000036
表示第j个sigma点对应协方差阵的权值，
Figure FDA0003488197360000037
表示矩阵P方根的第j列，λ＝α2(n+κ)-n用来降低总的预测误差，α控制采样点的分布状态，κ的选取应保证(n+λ)P是半正定矩阵，一般取值0，β≥0合并方程中高阶项动差；
3-2)、求状态一步预测及其对应的协方差阵：
Figure FDA0003488197360000038
Figure FDA0003488197360000039
Figure FDA00034881973600000310
其中：
Figure FDA00034881973600000311
是由于过程噪声统计量未知而根据经验选定的名义过程噪声协方差；
Figure FDA00034881973600000312
表示第j个sigma点的权值，
Figure FDA00034881973600000313
表示第j个sigma点对应协方差阵的权值。
5.根据权利要求1所述的一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：所述步骤4)具体包括以下步骤：
4-1)、选择IW分布作为一步预测协方差的先验分布：
Figure FDA00034881973600000314
其中，IW分布的自由度参数
Figure FDA00034881973600000315
逆尺度参数
Figure FDA00034881973600000316
通过如下方法得到：
对于一个IW分布：A～IW(A；t,T)，其期望可以写为：E[A-1]＝(t-n-1)Τ-1，t≥n+1，其中n是t的维度；因此，步骤3-2中的状态一步预测协方差Pk|k-1,l还可以表示为：
Figure FDA00034881973600000317
令：
Figure FDA0003488197360000041
则：
Figure FDA0003488197360000042
其中，nx是状态量x的维度，τl≥0是调谐参数，它的选取根据具体的情况而定；
4-2)、使用student’t分布作为量测的分布：
p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)
St(zk,l；hl(xk),Rk,l,vl)表示均值为hl(xk)，尺度矩阵为Rk,l，自由度参数为vl的student’t分布；各滤波器对目标状态估计独立进行，在各滤波器中，假设：p(zk,l|xk,l)＝St(zk,l；hl(xk,l),Rk,l,vl)≈p(zk,l|xk)＝St(zk,l；hl(xk),Rk,l,vl)，由于student’t分布的概率密度函数的闭合解难以求得，引入一个辅助随机变量λk,l，量测的概率密度函数可以写为如下积分形式：
Figure FDA0003488197360000043
其中，
Figure FDA0003488197360000044
表示形状参数为
Figure FDA0003488197360000045
和逆尺度参数为
Figure FDA0003488197360000046
的伽马分布；根据上一个式子，量测的概率密度函数最终可以表示为如下的分层高斯形式：
p(zk,l|xk,l,λk,l)＝N(zk,l；hl(xk,l),Rk,l/λk,l)
Figure FDA0003488197360000047
6.根据权利要求1所述的一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：所述步骤6)具体包括以下步骤：
6-1)、固定k时刻第i-1次迭代状态估计的近似后验概率分布q(i-1)(xk,l)，将第i次迭代的状态一步预测协方差Pk|k-1,l的近似后验概率分布更新为
Figure FDA0003488197360000048
其中：自由度参数
Figure FDA0003488197360000049
逆尺度参数
Figure FDA00034881973600000410
k时刻第i次迭代一步预测协方差的逆矩阵期望为
Figure FDA00034881973600000411
第i次迭代一步预测协方差可以表示为：
Figure FDA00034881973600000412
6-2)、固定k时刻第i-1次迭代状态估计的近似后验概率分布q(i-1)(xk,l)，将第i次迭代的辅助变量λk,l的近似后验概率分布更新为
Figure FDA00034881973600000413
其中：形状参数
Figure FDA0003488197360000051
逆尺度参数
Figure FDA0003488197360000052
Figure FDA0003488197360000053
Figure FDA0003488197360000054
是按照步骤3-1的方法对
Figure FDA0003488197360000055
进行无迹变换得到的第j个sigma点，
Figure FDA0003488197360000056
是对应的权值；
辅助变量期望
Figure FDA0003488197360000057
第i次迭代中辅助变量值表示为：
Figure FDA0003488197360000058
第i次迭代修正量测噪声协方差可以表示为：
Figure FDA0003488197360000059
6-3)、固定k时刻第i次迭代状态一步预测协方差Pk|k-1,l的近似后验概率分布q(i)(Pk|k-1,l)、第i次迭代辅助变量λk,l的近似后验概率分布q(i)(λk,l)，将第i次迭代的对目标状态估计xk,l的近似后验概率分布更新为
Figure FDA00034881973600000510
其中：
Figure FDA00034881973600000511
分别表示第l个滤波器在k时刻第i次迭代中的目标状态估计及其对应的协方差；
按照步骤3-1的方法对状态一步预测
Figure FDA00034881973600000512
和第i次迭代得到的状态一步预测协方差
Figure FDA00034881973600000513
进行无迹变换重新获得状态一步预测sigma点以及量测预测
Figure FDA00034881973600000514
Figure FDA00034881973600000515
Figure FDA00034881973600000516
其中：n表示
Figure FDA00034881973600000517
的维度，
Figure FDA00034881973600000518
表示通过无迹变换产生的第j个sigma点，
Figure FDA00034881973600000519
表示矩阵P方根的第j列，
Figure FDA00034881973600000520
表示对第j个sigma点的量测预测，
Figure FDA00034881973600000521
表示第j个sigma点对应的权值，
Figure FDA00034881973600000522
表示第l个滤波器第k时刻的量测预测；
在UKF框架下求得
Figure FDA00034881973600000523
Figure FDA00034881973600000524
Figure FDA00034881973600000525
Figure FDA0003488197360000061
Figure FDA0003488197360000062
Figure FDA0003488197360000063
7.根据权利要求1所述的一种基于变分贝叶斯的分布协同非线性系统状态估计方法，其特征在于：所述步骤8)具体包括以下步骤：
8-1)、k时刻局部滤波器状态估计的加权融合：
全局状态估计协方差：
Figure FDA0003488197360000064
全局状态估计：
Figure FDA0003488197360000065
8-2)、将全局估计及其协方差按照一定的准则反馈给各局部滤波器：
状态估计反馈：
Figure FDA0003488197360000066
状态估计协方差反馈：
Figure FDA0003488197360000067
其中，αk,l为反馈权重系数，它们随着各个局部滤波器协方差的变化而变化，αk,l满足以下条件：
αk,1+αk,2+...+αk,l＝1
Figure FDA0003488197360000068
其中：||·||F表示表示Frobenius范数，即对于任意矩阵A：
Figure FDA0003488197360000069
CN202210088496.7A 2022-01-25 2022-01-25 基于变分贝叶斯的分布协同非线性系统状态估计方法 Active CN114567288B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210088496.7A CN114567288B (zh) 2022-01-25 2022-01-25 基于变分贝叶斯的分布协同非线性系统状态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210088496.7A CN114567288B (zh) 2022-01-25 2022-01-25 基于变分贝叶斯的分布协同非线性系统状态估计方法

Publications (2)

Publication Number Publication Date
CN114567288A true CN114567288A (zh) 2022-05-31
CN114567288B CN114567288B (zh) 2024-04-26

Family

ID=81713773

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210088496.7A Active CN114567288B (zh) 2022-01-25 2022-01-25 基于变分贝叶斯的分布协同非线性系统状态估计方法

Country Status (1)

Country Link
CN (1) CN114567288B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117353705A (zh) * 2023-10-08 2024-01-05 哈尔滨工业大学 一种带有未知延迟概率的一步时延跟踪滤波方法及系统
CN117853827A (zh) * 2024-03-07 2024-04-09 安徽省大气探测技术保障中心 大气温室气体监测用采样泵工作状态运行监测系统及方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104112079A (zh) * 2014-07-29 2014-10-22 洛阳理工学院 一种模糊自适应变分贝叶斯无迹卡尔曼滤波方法
CN108599737A (zh) * 2018-04-10 2018-09-28 西北工业大学 一种变分贝叶斯的非线性卡尔曼滤波器的设计方法
CN108763167A (zh) * 2018-05-07 2018-11-06 西北工业大学 一种变分贝叶斯的自适应滤波方法
CN109508445A (zh) * 2019-01-14 2019-03-22 哈尔滨工程大学 一种带有色量测噪声和变分贝叶斯自适应卡尔曼滤波的目标跟踪方法
CN111178385A (zh) * 2019-12-02 2020-05-19 江苏大学 一种鲁棒在线多传感器融合的目标跟踪方法
WO2021237929A1 (zh) * 2020-05-27 2021-12-02 江南大学 基于贝叶斯学习的执行器故障估计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104112079A (zh) * 2014-07-29 2014-10-22 洛阳理工学院 一种模糊自适应变分贝叶斯无迹卡尔曼滤波方法
CN108599737A (zh) * 2018-04-10 2018-09-28 西北工业大学 一种变分贝叶斯的非线性卡尔曼滤波器的设计方法
CN108763167A (zh) * 2018-05-07 2018-11-06 西北工业大学 一种变分贝叶斯的自适应滤波方法
CN109508445A (zh) * 2019-01-14 2019-03-22 哈尔滨工程大学 一种带有色量测噪声和变分贝叶斯自适应卡尔曼滤波的目标跟踪方法
CN111178385A (zh) * 2019-12-02 2020-05-19 江苏大学 一种鲁棒在线多传感器融合的目标跟踪方法
WO2021237929A1 (zh) * 2020-05-27 2021-12-02 江南大学 基于贝叶斯学习的执行器故障估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YULONG HUANG等: "Robust Kalman Filters Based on Gaussian Scale Mixture Distributions With Application to Target Tracking", IEEE TRANS. SYSTEMS, MAN, AND CYBERNETICS: SYSTEMS, 15 December 2017 (2017-12-15) *
王荣;温阳;徐晓龙;: "一种新的基于闪烁噪声的扩展目标跟踪方法", 商洛学院学报, no. 02, 20 April 2016 (2016-04-20) *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117353705A (zh) * 2023-10-08 2024-01-05 哈尔滨工业大学 一种带有未知延迟概率的一步时延跟踪滤波方法及系统
CN117853827A (zh) * 2024-03-07 2024-04-09 安徽省大气探测技术保障中心 大气温室气体监测用采样泵工作状态运行监测系统及方法
CN117853827B (zh) * 2024-03-07 2024-05-14 安徽省大气探测技术保障中心 大气温室气体监测用采样泵工作状态运行监测系统及方法

Also Published As

Publication number Publication date
CN114567288B (zh) 2024-04-26

Similar Documents

Publication Publication Date Title
Wang et al. Zonotopic set-membership state estimation for discrete-time descriptor LPV systems
CN110647042B (zh) 一种基于数据驱动的机器人鲁棒学习预测控制方法
EP1015943B1 (en) A method for real-time nonlinear system state estimation and control
CN108132599B (zh) 一种基于迭代反馈整定的ude控制系统设计方法
CN109827579B (zh) 一种组合定位中滤波模型实时校正的方法和系统
CN110609476B (zh) 一种基于高斯过程模型的多变量非线性动态系统模型预测控制方法
CN110826021B (zh) 一种非线性工业过程鲁棒辨识和输出估计方法
CN116088307B (zh) 基于误差触发自适应稀疏辨识的多工况工业过程预测控制方法、装置、设备及介质
Zhao et al. Robust FIR state estimation of dynamic processes corrupted by outliers
CN114567288A (zh) 基于变分贝叶斯的分布协同非线性系统状态估计方法
Jeong et al. Adaptive event-triggered tracking using nonlinear disturbance observer of arbitrarily switched uncertain nonlinear systems in pure-feedback form
CN114815619A (zh) 一种量测随机丢失下基于卡尔曼滤波的机器人跟踪方法
CN113452349A (zh) 一种基于贝叶斯序贯重要性积分的卡尔曼滤波方法
Malik et al. State and output feedback local control schemes for nonlinear discrete-time 2-D Roesser systems under saturation, quantization and slope restricted input
Wolf et al. Rigorous solution vs. fast update: Acceptable computational delay in NMPC
Zamani et al. Minimum-energy filtering on the unit circle
Bonzanini et al. On the stability properties of perception-aware chance-constrained mpc in uncertain environments
CN110826184B (zh) 一种在时变时滞下narx模型结构和参数的变分贝叶斯辨识方法
Kim et al. A single-pass noise covariance estimation algorithm in adaptive Kalman filtering for non-stationary systems
Kim et al. A single-pass noise covariance estimation algorithm in multiple-model adaptive kalman filtering for non-stationary systems
Jiang et al. Fast and smooth composite local learning-based adaptive control
Chalupa Predictive control using self tuning model predictive controllers library
Herzallah et al. Robust control of nonlinear stochastic systems by modelling conditional distributions of control signals
AU2021102341A4 (en) A model order reduction method for optimal model approximation of linear-time invariant systems
Herzallah et al. Distribution modeling of nonlinear inverse controllers under a Bayesian framework

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