卡尔曼滤波的五大公式及python代码示例_卡尔曼滤波器五个公式-程序员宅基地

技术标签: python  计算机视觉|目标跟踪  开发语言  

一、系统的状态方程

状态方程是根据上一时刻的状态和控制变量来推测此刻的状态
在这里插入图片描述

  • x k x_{k} xk 是状态分量的 n n n维矢量
  • A A A 是  n ∗ n n * n nn 的状态转移矩阵,也就是对目标状态转换的猜想模型,是已知的
  • u k − 1 u_{k-1} uk1 是新的,让系统可以接受外部控制
  • B B B 是  n ∗ c n * c nc 矩阵,将输入转换为状态的矩阵
  • w k − 1 w_{k-1} wk1 是预测过程的噪声,对应 x k x_{k} xk中每个分量的噪声,期望为0,协方差为 Q Q Q的高斯白噪声

二、观测方程

在这里插入图片描述

  • H H H 是  m ∗ n m * n mn矩阵,是状态变量(观测)的转换矩阵,表示将状态和观测连接起来的关系,卡尔曼滤波里为线性关系,它负责将 m 维的测量值转换到 n 维,使之符合状态变量的数学形式,是滤波的前提条件之一。
  • v k v_{k} vk 观测噪声,服从高斯分布 N ( 0 , R ) N(0, R) N(0,R) , R R R 即下文测量噪声

三、五大公式

以下摘自
卡尔曼滤波五个公式各个参数的意义

  • 信息过程足够精确的模型,是由白噪声所激发的线性(也可以是时变)的动态系统
  • 每次测量的信号都包含这附加的白噪声分量
    当满足以上假设时,可以应用卡尔曼滤波算法

卡尔曼滤波算法有两个基本假设:

卡尔曼滤波算法分为两步:预测和更新

  • 预测 : 根据上一时刻( k - 1 k - 1 k1 时刻) 的后验估计值来估计当前时刻( k k k 时刻) 的状态,得到 k k k 时刻的先验估计值;
  • 更新: 使用当前时刻的测量值来更正预测阶段估计值,得到当前时刻的后验估计值。
    卡尔曼滤波器可以分为时间更新方程和测量更新方程。时间更新方程(即预测阶段)根据前一时刻的状态估计值推算当前时刻的状态变量先验估计值和误差协方差先验估计值; 测量更新方程(即更新阶段)负责将先验估计和新的测量变量结合起来构造改进的后验估计。时间更新方程和测量更新方程也被称为预测方程和校正方程。因此卡尔曼算法是一个递归的预测—校正方法。

在这里插入图片描述

  • 1,: 分别表示 k - 1 时刻和 k 时刻的后验状态估计值,是滤波的结果之一,即更新后的结果,也叫最优估计(估计的状态,根据理论,我们不可能知道每时刻状态的确切结果所以叫估计)。

  • 2,: k 时刻的先验状态估计值,是滤波的中间计算结果,即根据上一时刻(k-1时刻)的最优估计预测的k时刻的结果,是预测方程的结果。

  • 3,分别表示 k - 1 时刻和 k 时刻的后验估计协方差(即的协方差,表示状态的不确定度),是滤波的结果之一。

  • 4,k 时刻的先验估计协方差(的协方差),是滤波的中间计算结果。

  • 5,是状态变量到测量(观测)的转换矩阵,表示将状态和观测连接起来的关系,卡尔曼滤波里为线性关系,它负责将 m 维的测量值转换到 n 维,使之符合状态变量的数学形式,是滤波的前提条件之一。

  • 6,测量值(观测值),是滤波的输入。

  • 7,滤波增益矩阵,是滤波的中间计算结果,卡尔曼增益,或卡尔曼系数。

  • 8,状态转移矩阵,实际上是对目标状态转换的一种猜想模型。例如在机动目标跟踪中, 状态转移矩阵常常用来对目标的运动建模,其模型可能为匀速直线运动或者匀加速运动。当状态转移矩阵不符合目标的状态转换模型时,滤波会很快发散。

  • 9,Q:过程激励噪声协方差(系统过程的协方差)。该参数被用来表示状态转换矩阵与实际过程之间的误差。因为我们无法直接观测到过程信号, 所以 Q 的取值是很难确定的。是卡尔曼滤波器用于估计离散时间过程的状态变量,也叫预测模型本身带来的噪声。状态转移协方差矩阵
     

  • 10:R: 测量噪声协方差。滤波器实际实现时,测量噪声协方差 R一般可以观测得到,是滤波器的已知条件。

  • 11,B:是将输入转换为状态的矩阵

四、Python代码

这里真实值为 x=-0.377,且假设A=1, H=1
观测值存在噪声,那么如何估计出实际的值呢?
这里给出两种方案,一种是北卡大学开源的,直接通过公式计算的结果

# -*- coding=utf-8 -*-  
  # Kalman filter example demo in Python  
     
   # A Python implementation of the example given in pages 11-15 of "An  
   # Introduction to the Kalman Filter" by Greg Welch and Gary Bishop,  
   # University of North Carolina at Chapel Hill, Department of Computer  
   # Science, TR 95-041,  
   # http://www.cs.unc.edu/~welch/kalman/kalmanIntro.html  
     
   # by Andrew D. Straw  
   #coding:utf-8  
   import numpy  
   import pylab  
     
   #这里是假设A=1,H=1的情况  
     
   # 参数初始化  
   n_iter = 50  
   sz = (n_iter,) # size of array  
   x = -0.37727 # 真实值  
   z = numpy.random.normal(x,0.1,size=sz) # 观测值 ,观测时存在噪声
     
   Q = 1e-5 # process variance  
     
   # 分配数组空间  
   xhat=numpy.zeros(sz)      # x 滤波估计值  
   P=numpy.zeros(sz)         # 滤波估计协方差矩阵  
   xhatminus=numpy.zeros(sz) #  x 估计值  
   Pminus=numpy.zeros(sz)    # 估计协方差矩阵  
   K=numpy.zeros(sz)         # 卡尔曼增益  
     
   R = 0.1**2 # estimate of measurement variance, change to see effect  
     
   # intial guesses  
   xhat[0] = 0.0  
   P[0] = 1.0  
     
   for k in range(1,n_iter):  
       # 预测  
       xhatminus[k] = xhat[k-1]  #X(k|k-1) = AX(k-1|k-1) + BU(k) + W(k),A=1,BU(k) = 0  
       Pminus[k] = P[k-1]+Q      #P(k|k-1) = AP(k-1|k-1)A' + Q(k) ,A=1  
     
       # 更新  
       K[k] = Pminus[k]/( Pminus[k]+R ) #Kg(k)=P(k|k-1)H'/[HP(k|k-1)H' + R],H=1  
       xhat[k] = xhatminus[k]+K[k]*(z[k]-xhatminus[k]) #X(k|k) = X(k|k-1) + Kg(k)[Z(k) - HX(k|k-1)], H=1  
       P[k] = (1-K[k])*Pminus[k] #P(k|k) = (1 - Kg(k)H)P(k|k-1), H=1  
     
   pylab.figure()  
   pylab.plot(z,'k+',label='noisy measurements')     #观测值  
   pylab.plot(xhat,'b-',label='a posteri estimate')  #滤波估计值  
   pylab.axhline(x,color='g',label='truth value')    #真实值  
   pylab.legend()  
   pylab.xlabel('Iteration')  
   pylab.ylabel('Voltage')  
     
   pylab.figure()  
   valid_iter = range(1,n_iter) # Pminus not valid at step 0  
   pylab.plot(valid_iter,Pminus[valid_iter],label='a priori error estimate')  
   pylab.xlabel('Iteration')  
   pylab.ylabel('$(Voltage)^2$')  
   pylab.setp(pylab.gca(),'ylim',[0,.01])  
   pylab.show()  

另外一种是调用from filterpy.kalman里的卡尔曼滤波函数

from filterpy.kalman import KalmanFilter
import numpy as np
np.random.seed(0)
kf = KalmanFilter(dim_x=1, dim_z=1)
kf.F = np.array([1])
kf.H = np.array([1])
kf.R = np.array([0.1**2])
kf.P = np.array([1.0])
kf.Q = 1e-5 
xhat[0] = 0.0  
P[0] = 1.0 
for k in range(1,n_iter):  
    kf.predict()
    xhat[k] = kf.x
    kf.update(z[k], 0.1**2, np.array([1]))
    
    
    

    
pylab.figure()  
pylab.plot(z,'k+',label='noisy measurements')     #观测值  
pylab.plot(xhat,'b-',label='a posteri estimate')  #滤波估计值  
pylab.axhline(x,color='g',label='truth value')    #真实值  
pylab.legend()  
pylab.xlabel('Iteration')  
pylab.ylabel('Voltage')  

pylab.figure()  
valid_iter = range(1,n_iter) # Pminus not valid at step 0  
pylab.plot(valid_iter,Pminus[valid_iter],label='a priori error estimate')  
pylab.xlabel('Iteration')  
pylab.ylabel('$(Voltage)^2$')  
pylab.setp(pylab.gca(),'ylim',[0,.01])  
pylab.show()  

python-opencv中的卡尔曼滤波函数

kalman = cv2.KalmanFilter(1, 1)
kalman.transitionMatrix = np.array([[1]], np.float32)  # 转移矩阵 A
kalman.measurementMatrix = np.array([[1]], np.float32)  # 测量矩阵    H
kalman.measurementNoiseCov = np.array([[1]], np.float32) * 0.01 # 测量噪声        R
kalman.processNoiseCov = np.array([[1]], np.float32) * 1e-5  # 过程噪声 Q
kalman.errorCovPost = np.array([[1.0]], np.float32)  # 最小均方误差 P

xhat = np.zeros(sz)  # x 滤波估计值 
kalman.statePost = np.array([xhat[0]], np.float32)
for k in range(1, n_iter):
#     print(np.array([z[k]], np.float32))
    mes = np.reshape(np.array([z[k]], np.float32), (1, 1))
# #     print(mes.shape)
    xhat[k] = kalman.predict()
    kalman.correct(np.array(mes, np.float32))



pylab.figure()
pylab.plot(z, 'k+', label='noisy measurements')  # 观测值
pylab.plot(xhat, 'b-', label='a posteri estimate')  # 滤波估计值
pylab.axhline(x, color='g', label='truth value')  # 真实值
pylab.legend()
pylab.xlabel('Iteration')
pylab.ylabel('Voltage')
pylab.show() 

三者都能得到同一结果
在这里插入图片描述

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/hjxu2016/article/details/114102906

智能推荐

Linux 磁盘管理及磁盘分区(主分区、扩展分区、逻辑分区)_linux分区-程序员宅基地

文章浏览阅读9.7k次,点赞26次,收藏104次。Linux 磁盘管理及磁盘分区(主分区、扩展分区、逻辑分区)MBR xfsswapFAT16 FAT32EXT4JFS_linux分区

java关键字abstract(抽象)详解_abstractjava-程序员宅基地

文章浏览阅读10w+次,点赞227次,收藏778次。说到抽象,想必大家第一个联想到的就是抽象画,基本上看不懂画的是啥,只能说它在表达某个概念。 在java中我们中abstract关键字来表达抽象。举个例子: 我们说车子都可以跑(run)。但有几个轮子,怎么跑,对于不同的车有不同的结果。自行车需要人踩着跑,汽车发动机推动跑等等,那么我们可以车表达为抽象类。/** * 车子类 */public abstract class Car {_abstractjava

使用el-table选择进行保存,实现分页多选,设置选中,编辑原数据默认值不变_el-table 多选框默认选中 样式不更新-程序员宅基地

文章浏览阅读1.7k次。el-table分页选择,默认选中设置,二次选择原商品编辑的数据不变_el-table 多选框默认选中 样式不更新

信息安全产品认证_网络关键设备和网络安全专用产品安全认证证书-程序员宅基地

文章浏览阅读5.6k次。目前市场上通行的信息安全产品认证有三种,分别是:1.《网络关键设备和网络安全专用产品安全认证证书》2.《中国国家信息安全产品认证证书》3.《IT产品信息安全认证证书》其中《中国国家信息安全产品认证证书》、《网络关键设备和网络安全专用产品安全认证证书》为强制性认证证书,在产品功能和性能满足条件的情况下,两张证书可以同一产品同时认证同时发证。《IT产品信息安全认证证书》为自愿性认证证书。_网络关键设备和网络安全专用产品安全认证证书

从零构建一个微信小程序_从零搭建一个微信小程序-程序员宅基地

文章浏览阅读1.3k次。一、微信小程序介绍微信小程序,简称小程序,英文名:Mini Program,是一种不需要下载安装即可使用的应用,它实现了应用“触手可及”的梦想,用户扫一扫或搜一下即可打开应用。1.1 为什么是微信小程序?微信有海量用户,而且粘性很高,在微信里开发产品更容易触达用户。推广 app 或公众号的成本太高。开发适配成本低。容易小规模试错,然后快速迭代。跨平台。1.2 还有其他的小程序 不容忽视支付宝小程序百度小程序QQ 小程序今日头条 + 抖音小程序1.3 优秀的第三方小程序_从零搭建一个微信小程序

前端开发:CSS的“*”“#”“.”符号的对比使用_css中".","#","*"的区别和作用-程序员宅基地

文章浏览阅读4.9k次,点赞2次,收藏7次。在前端开发过程中,CSS是前端的三驾马车之一,关于CSS的使用也是非常重要的环节,可能有人觉得在前端学习或开发中最容易的就是关于CSS的,但是要想精通使用CSS也不是一件容易的事情,也就是看着简单要想完美使用难,所以CSS的使用也是不容忽视的。本篇博文来分享一个关于CSS的符号相关内容,关于它的使用也是比较常用的,常用的有:“*”、“#”、“.”等符号,所以在前端项目开发中经常会用到该功能,总结一下,方便查阅使用。.................._css中".","#","*"的区别和作用

随便推点

Unity shader 实现图片带圆角和边线border_unity图片圆角-程序员宅基地

文章浏览阅读4.7k次,点赞6次,收藏19次。unity shader实现圆角+ border_unity图片圆角

Python数据处理040:数据分析之Excel文件_python分析excel数据-程序员宅基地

文章浏览阅读2.9k次。Python数据分析之Excel文件_python分析excel数据

11. python入门复习教程之命名空间与作用域,再谈异常,标准库,持久化与序列化_websocket 命名空间的作用-程序员宅基地

文章浏览阅读1.9w次,点赞10次,收藏14次。如果你学习 python 一段时间之后,会经常碰到命名空间,作用域 这两个名词,接下来咱们就从这两个名词开始入手,展开这篇博客。命名空间与作用域命名空间是指对象所处的位置,命名二字一般指的是变量名,举一个例子,例如在某个函数内,如果使用了 a = 100 这行代码,遵循的规则是在这个命名空间内,如果之前没有变量 a,那此时定义 a 并且对 a 赋值,如果之前存在 a,那对之前的赋值进行覆盖,这个函数内部,就是变量 a 的命名空间。作用域指的是某个对象的有效作用范围,这个作用范围又由命名空间决定,所以在_websocket 命名空间的作用

转:Android原生PDF功能实现:PDF阅读、PDF页面跳转、PDF手势伸缩、PDF目录树、PDF预览缩略图_android pdf 目录-程序员宅基地

文章浏览阅读511次。基于Google的android-pdfview拓展的pdf阅读功能:(先把地址保存下来,就怕大佬哪天心情不好给删了)简书大佬博客:https://www.jianshu.com/p/eec432fa89afDemo下载地址(github + 百度网盘):https://github.com/qxcwanxss/AndroidPdfViewerDemohttps://pan.baidu.c..._android pdf 目录

关于app inventor跨界面连接蓝牙并且一键SOS拨打电话求救(三)_app inventor 电话拨号器-程序员宅基地

文章浏览阅读2.5k次,点赞3次,收藏18次。之前写的两个博客都是在做一个项目中遇到的问题,因为自己没有在网上找到太多资料,所以记录了下来,第三个博客是把前两个博客连接了起来。有不懂得可以看看前面两个,或者在下面留言。Arduino和蓝牙程序和第一个博客的程序一样,只是改了一下延时const int buttonPin = 2; // 连接按键的引脚 // 变量定义int buttonState = 0; ..._app inventor 电话拨号器

超好用的几个开发工具使用_sqlitestudio自动补全-程序员宅基地

文章浏览阅读8.3k次。1. Makefile 作用2. Makefile的语法规则3. Makefile使用_sqlitestudio自动补全

推荐文章

热门文章

相关标签