卡尔曼滤波的五大公式及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

智能推荐

互联网广告及产品变现认知分析整理_大厂商业化变现部门会给佣金吗-程序员宅基地

文章浏览阅读1.4w次,点赞5次,收藏6次。无论是互联网企业还是传统的其他企业情况,变现模式比如说用户直接付费、佣金分成、内容付费、电子商务、金融运作、增值服务等,哪怕刚开始选择的是其他变现模式,最终也还是免不了要增加广告产品进行变现,也就是说广告才是互联网产品最常见的变现模式。本文对广告产品进行基本分析和变现模式分析_大厂商业化变现部门会给佣金吗

使用Doxygen把代码文档化_doxygen后显示无法被文档化-程序员宅基地

文章浏览阅读96次。https://www.fgba.net/sitemap.xml_doxygen后显示无法被文档化

Cocos2d-x在win7下的android交叉编译环境-程序员宅基地

文章浏览阅读763次,点赞8次,收藏17次。很多人在刚接触这个行业的时候或者是在遇到瓶颈期的时候,总会遇到一些问题,比如学了一段时间感觉没有方向感,不知道该从哪里入手去学习,对此我整理了一些资料,需要的可以免费分享给大家我的【Github】会分享一些关于Android进阶方面的知识,也会分享一下最新的面试题~如果你熟练掌握GitHub中列出的知识点,相信将会大大增加你通过前两轮技术面试的几率!这些内容都供大家参考,互相学习。

轮番面试,阿里P8总结出了java程序设计,Javaweb框架面试题-程序员宅基地

文章浏览阅读256次,点赞5次,收藏10次。由于文案过于长,在此就不一一介绍了,这份Java后端架构进阶笔记内容包括:Java集合,JVM、Java并发、微服务、SpringNetty与 RPC 、网络、日志 、Zookeeper 、Kafka 、RabbitMQ 、Hbase 、MongoDB、Cassandra 、Java基础、负载均衡、数据库、一致性算法、Java算法、数据结构、分布式缓存等等知识详解。本知识体系适合于所有Java程序员学习,关于以上目录中的知识点都有详细的讲解及介绍,掌握该知识点的所有内容对你会有一个质的提升,

Android~apk的混淆和加固(1),你有过迷茫吗-程序员宅基地

文章浏览阅读579次,点赞29次,收藏11次。针对apk,加固是多维度的安全防护方案,包括反破解、反逆向、防篡改等,可以防止应用被各类常见破解工具逆向,安全性要远大于单纯的代码混淆。操作的对象是项目打包成的apk文件。为什么我们需要混淆?因为java字节码特性很容易反编译。对于加固,上架应用市场一般提供相关文档指导我们进行apk的渠道打包发布,这里不做展开我们先大概知道加固的一些原理。修改gradlerelease {参考需要保留的类及方法,确定项目中哪些不能混淆的类参考混淆模板,编写我们的混淆文件。

RAS使用拨号网络拨号的类_mfc ras拨号不用指定域名吗-程序员宅基地

文章浏览阅读4.7k次。RAS Socket 拨号网络_mfc ras拨号不用指定域名吗

随便推点

PostgreSQL的存储过程及基本使用_postgre 调用存储过程-程序员宅基地

文章浏览阅读5k次。PostgreSQL的存储过程及基本使用 一、存储过程的结构二、变量使用1.变量类型2.record变量3.赋值 三、基本流程语句1.分支选择2.循环 四、查询并返回多条记录五、其它 一、存储过程的结构 一个求长方形面积的存储过程。(当然,这个存储过程在数据库中并没有什么实用意义,这..._postgre 调用存储过程

彻底弄清楚C++ static 静态成员与静态成员函数的原理_static cpp 成员函数-程序员宅基地

文章浏览阅读634次。参考mooc魏英《C++程序设计》文字是魏老师的讲解,纯手打,ppt上没有。为什么会用到静态成员:现在大型应用程序都是由多个程序员所开发的,那么多个程序员就需要使用一个共同都能使用的数据来解决一些问题,采用静态数据成员解决这一问题。思考一下现在我们需要统计员工的总人数,能不能在这个员工类中增加一个成员专门用来存放总人数呢?这样做是不好的:1.每个对象都要增加一个这样的成员,对存储空间是一种浪费。对于公司来说,总人数只有一个值,那么每个对象都要增加一个这样的成员,浪费存储空间。2.使用不方便,当总人数_static cpp 成员函数

Qt上报错:undefined reference to xxx@openssl_1.0.0_undefined reference to symbol 'x509_get_subject_na-程序员宅基地

文章浏览阅读3.8k次,点赞2次,收藏6次。关于openssl版本不兼容问题产生原因解决办法原本pro配置文件:/*无效的配置*/#增加程序库文件路径LIBS += \ -lcurl -lcrypto -ljsoncpp直接上错误代码:/usr/lib/x86_64-linux-gnu/libcurl.so:对‘PKCS12_PBE_add@OPENSSL_1.0.0’未定义的引用/usr/lib/x86_64-linux-gnu/libcurl.so:对‘OCSP_basic_verify@OPENSSL__undefined reference to symbol 'x509_get_subject_name@@openssl_1_1_0

N叉树(N-ary Tree)-程序员宅基地

文章浏览阅读2.7k次。N叉树N 叉树的定义N 叉树的遍历N 叉树的经典递归解法"自顶向下"的解决方案"自底向上"的解决方案N 叉树的定义二叉树是一棵以根节点开始,每个节点含有不超过 2 个子节点的树。将这个定义扩展到 N 叉树 。 一棵以根节点开始,每个节点不超过 N 个子节点的树,称为 N叉树阅读参考:https://leetcode-cn.com/leetbook/read/n-ary-tree/x0wi57/N 叉树的遍历一棵二叉树可以按照前序、中序、后序或者层序来进行遍历。在这些遍历方法中,前序遍历、后序遍历和_n叉树

shell脚本编程常见系统变量和环境变量,shell脚本常用流程控制语句_熟悉shell常用的系统变量和环境变量,练习语句-程序员宅基地

文章浏览阅读6.1k次。shell编程常见系统变量命令功能$0当前脚本名称$?命令或程序执行完后的状态,返回0则表示执行成功$n当前脚本的第n个参数,n=1,2,…9$#当前脚本的参数个数(不包括程序本身)$*当前脚本的所有参数(不包括程序本身)$$程序本身的PID号shell编程常见环境变量命令功能PATH命令所示路径,以冒号为分割PWD显示当前所在路径USER打印当前用户名ID打印当前用户id信息TERM_熟悉shell常用的系统变量和环境变量,练习语句

PostgreSQL 获取拼音首字母的函数 (可以获取所有中文字符)经典原创分享-程序员宅基地

文章浏览阅读1.4k次。CREATE OR REPLACE FUNCTION "gis"."get_hzpycap"("v_str" varchar, "needhz" int4=0) RETURNS "pg_catalog"."varchar" AS $BODY$DECLARE pos INT4; hzlen INT4; hz ..._pg 中文简拼函数

推荐文章

热门文章

相关标签