草庐IT

滤波笔记一:卡尔曼滤波(Kalman Filtering)详解

scoutee 2023-04-24 原文

本笔记是总结了B站DR_CAN的卡尔曼滤波器的课程,他的B站主页为:DR_CAN的个人空间_哔哩哔哩_bilibili

PS:虽然我不是学自控的,但是老师真的讲的很好! 

目录

Lesson1 递归算法

Lesson2 数学基础_数据融合_协方差矩阵_状态空间方程

Lesson3 卡尔曼增益的详细推导

Lesson4 误差的协方差矩阵Pe的数学推导

 Lesson5 直观理解卡尔曼滤波以及一个实例

当计算误差Wk大于测量误差Vk时

当计算误差Wk小于测量误差Vk时

本例的python代码

突然想到一个问题:如何确定卡尔曼滤波要迭代多少次呢?

总结一下

1.算法迭代的五个步骤

2.算法的python代码实现


Lesson1 递归算法

 

 

Lesson2 数学基础_数据融合_协方差矩阵_状态空间方程

 

 

Lesson3 卡尔曼增益的详细推导

 

 

 

Lesson4 误差的协方差矩阵Pe的数学推导

 

 

 Lesson5 直观理解卡尔曼滤波以及一个实例

 

 

 

 下面具体看一下,之前反复提到过:如果模型计算误差Wk小,最终的估计值就更偏向计算值;

                                                           如果测量误差Vk小,最终的估计值就更偏向于测量值。

而在这个示例中,Wk/Vk的偏差是以其协方差矩阵来反映的(主对角线是方差)。

当计算误差Wk大于测量误差Vk时

Wk的协方差矩阵为Q,Vk的协方差矩阵为R:

 结果图:

结果分析:

真实的实际速度是蓝色曲线,最终的估计速度即为后验估计速度,是黄色曲线。对比它们之间的偏差能够看出估计值与实际值之间的误差,从而判断算法的准确度。

在最优化方法中我们知道:最优估计有不同的准则,比如:最小二乘估计、最小方差估计、极大后验估计等等。具体内容不赘述。

我们要知道,如果没有不确定性(即Wk和Vk),那么估计值就是实际值(精准估计)。

卡尔曼滤波中采用的就是使得误差的方差最小为最优估计准则因为如果后验估计值和真实值越接近,那么误差ek的变化就很小,即误差ek的方差很小。

进一步推导,考虑到误差ek会有很多不同的分量(因为状态量不同,比如说此例子中就是有状态量X1表示位置,状态量X2表示速度,那么它们就分别有误差e1和e2)。要使得总误差方差最小,那么误差各个分量的方差之和加起来就要最小。而“误差各分量的方差之和”正好是误差的协方差矩阵的主对角线之和——迹。

故此我们引入了Wk的协方差矩阵为Q,Vk的协方差矩阵为R。然后其方差越大时,说明误差越大,即越不可以相信。所以此处计算误差较大,可以看到先验估计速度(灰色)偏离实际速度(蓝色)的程度要大于测量速度(红色)偏离实际速度(蓝色的程度)。

所以最后的估计值——后验估计速度(黄色)曲线也是更为接近测量速度(红色)曲线。

一种通俗的理解方式就是:建模计算值和测量值都是不准确的,两者的不准确程度分别以计算误差的方差测量误差的方差来衡量,方差越大越不可以相信。在两个不准确的值的基础上尽量准确估计,就是谁方差越小,越相信谁,越靠近谁。

当计算误差Wk小于测量误差Vk时

 结果分析:

由于此时测量误差的方差较大,导致测量值很不可信,其变化的程度可以看到也很离谱。但是由于后验估计值(黄色)更为依赖模型计算值(灰色),所以后验估计值也没离实际值(蓝色)太远。

而这,正是卡尔曼滤波的作用:在不准确之中得到最准确的估计值。

本例的python代码

源代码来源:B站用户东爱北的GitHub​​​​​​https://github.com/liuchangji/2D-Kalman-Filter-Example_Dr_CAN_in_python

我对其中的源码做了注释,以及对一个小错误进行了修改(产生符合高斯分布的变量时,scale输入的应该是标准差,而协方差矩阵里面主对角线上面是方差,所以要开根号,要注意开完根号要保证其类型仍为np.float) 



import numpy as np
import math
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号


#定义一个产生符合高斯分布的函数,均值为loc=0.0,标准差为scale=sigma,输出的大小为size
def gaussian_distribution_generator(sigma):
    return np.random.normal(loc=0.0,scale=sigma,size=None)

# 状态转移矩阵,上一时刻的状态转移到当前时刻
A = np.array([[1,1],[0,1]])

# 过程噪声w协方差矩阵Q,P(w)~N(0,Q),噪声来自真实世界中的不确定性
Q = np.array([[0.01,0],[0,0.01]])

# 测量噪声协方差矩阵R,P(v)~N(0,R),噪声来自测量过程的误差
R = np.array([[1,0],[0,1]])

# 传输矩阵/状态观测矩阵H
H = np.array([[1,0],[0,1]])

# 控制输入矩阵B
B = None

# 初始位置和速度
X0 = np.array([[0],[1]])

# 状态估计协方差矩阵P初始化
P =np.array([[1,0],[0,1]])

if __name__ == "__main__":
    #---------------------初始化-----------------------------
    #真实值初始化 这里还要再写一遍np.array是保证它的类型是数组array
    X_true = np.array(X0)
    #后验估计值Xk的初始化
    X_posterior = np.array(X0)
    #第k次误差的协方差矩阵的初始化
    P_posterior = np.array(P)

    #创建状态变量的真实值的矩阵 状态变量1:速度 状态变量2:位置
    speed_true = []
    position_true = []

    #创建测量值矩阵
    speed_measure = []
    position_measure = []

    #创建状态变量的先验估计值
    speed_prior_est = []
    position_prior_est = []

    #创建状态变量的后验估计值
    speed_posterior_est = []
    position_posterior_est = []

    #---------------------循环迭代-----------------------------
    #设置迭代次数为30次

    for i in range(30):
        #--------------------建模真实值-----------------------------

        # 生成过程噪声w w=[w1,w2].T(列向量)
        # Q[0,0]是过程噪声w的协方差矩阵的第一行第一列,即w1的方差,Q[1,1]是协方差矩阵的第二行第二列,即为w2的方差
        # python的np.random.normal(loc,scale,size)函数中scale输入的是标准差,所以要开方
        Q_sigma = np.array([[math.sqrt(Q[0,0]),Q[0,1]],[Q[1,0],math.sqrt(Q[1,1])]])
        w = np.array([[gaussian_distribution_generator(Q_sigma[0, 0])],
                      [gaussian_distribution_generator(Q_sigma[1, 1])]])
        # print('00',Q[0,0],'它的类型是',type(Q[0,0]))
        # print('开根号的00', Q_sigma[0, 0], '它的类型是', type(Q_sigma[0, 0]))
        # print('00的平方根',math.sqrt(Q[0,0]),"它的类型是",type(math.sqrt(Q[0,0])))
        # print('w[',i,']=',w)

        # 真实值X_true 得到当前时刻的状态;之前我一直在想它怎么完成从Xk-1到Xk的更新,实际上在代码里面直接做迭代就行了,这里是不涉及数组下标的!!!
        #dot函数用于矩阵乘法,对于二维数组,它计算的是矩阵乘积
        X_true = np.dot(A, X_true) + w

        # 速度的真实值是speed_true 使用append函数可以把每一次循环中产生的拼接在一起,形成一个新的数组speed_true
        speed_true.append(X_true[1,0])
        position_true.append(X_true[0,0])
        #print(speed_true)


        # --------------------生成观测值-----------------------------
        # 生成过程噪声
        R_sigma = np.array([[math.sqrt(R[0,0]),R[0,1]],[R[1,0],math.sqrt(R[1,1])]])
        v = np.array([[gaussian_distribution_generator(R_sigma[0,0])],[gaussian_distribution_generator(R_sigma[1,1])]])

        # 生成观测值Z_measure 取H为单位阵
        Z_measure = np.dot(H, X_true) + v
        speed_measure.append(Z_measure[1,0]),
        position_measure.append(Z_measure[0,0])

        # --------------------进行先验估计-----------------------------
        # 开始时间更新
        # step1:基于k-1时刻的后验估计值X_posterior建模预测k时刻的系统状态先验估计值X_prior
        # 此时模型控制输入U=0
        X_prior = np.dot(A, X_posterior)
        # 把k时刻先验预测值赋给两个状态分量的先验预测值 speed_prior_est = X_prior[1,0];position_prior_est=X_prior[0,0]
        # 再利用append函数把每次循环迭代后的分量值拼接成一个完整的数组
        speed_prior_est.append(X_prior[1,0])
        position_prior_est.append(X_prior[0,0])

        # step2:基于k-1时刻的误差ek-1的协方差矩阵P_posterior和过程噪声w的协方差矩阵Q 预测k时刻的误差的协方差矩阵的先验估计值 P_prior
        P_prior_1 = np.dot(A, P_posterior)
        P_prior = np.dot(P_prior_1, A.T) + Q

        # --------------------进行状态更新-----------------------------
        # step3:计算k时刻的卡尔曼增益K
        k1 = np.dot(P_prior, H.T)
        k2 = np.dot(H, k1) + R
        #k3 = np.dot(np.dot(H, P_prior), H.T) + R  k2和k3是两种写法,都可以
        K = np.dot(k1, np.linalg.inv(k2))

        # step4:利用卡尔曼增益K 进行校正更新状态,得到k时刻的后验状态估计值 X_posterior
        X_posterior_1 = Z_measure -np.dot(H, X_prior)
        X_posterior = X_prior + np.dot(K, X_posterior_1)
        # 把k时刻后验预测值赋给两个状态分量的后验预测值 speed_posterior_est = X_posterior[1,0];position_posterior_est = X_posterior[0,0]
        speed_posterior_est.append(X_posterior[1,0])
        position_posterior_est.append(X_posterior[0,0])

        # step5:更新k时刻的误差的协方差矩阵 为估计k+1时刻的最优值做准备
        P_posterior_1 = np.eye(2) - np.dot(K, H)
        P_posterior = np.dot(P_posterior_1, P_prior)

    # ---------------------再从step5回到step1 其实全程只要搞清先验后验 k的自增是隐藏在循环的过程中的 然后用分量speed和position的append去记录每一次循环的结果-----------------------------

    # 可视化显示 画出速度比较和位置比较
    if True:
        # 画出1行2列的多子图
        fig, axs = plt.subplots(1,2)
        #速度
        axs[0].plot(speed_true,"-",color="blue",label="速度真实值",linewidth="1")
        axs[0].plot(speed_measure,"-",color="grey",label="速度测量值",linewidth="1")
        axs[0].plot(speed_prior_est,"-",color="green",label="速度先验估计值",linewidth="1")
        axs[0].plot(speed_posterior_est,"-",color="red",label="速度后验估计值",linewidth="1")
        axs[0].set_title("speed")
        axs[0].set_xlabel('k')
        axs[0].legend(loc = 'upper left')


        #位置
        axs[1].plot(position_true,"-",color="blue",label="位置真实值",linewidth="1")
        axs[1].plot(position_measure,"-",color="grey",label="位置测量值",linewidth="1")
        axs[1].plot(position_prior_est,"-",color="green",label="位置先验估计值",linewidth="1")
        axs[1].plot(position_posterior_est,"-",color="red",label="位置后验估计值",linewidth="1")
        axs[1].set_title("position")
        axs[1].set_xlabel('k')
        axs[1].legend(loc = 'upper left')

        #     调整每个子图之间的距离
        plt.tight_layout()
        plt.figure(figsize=(60, 40))
        plt.show()

 结果图1(迭代30次):

图2(迭代60次):

 图1结果分析:

本次实例中,取了过程噪声的协方差矩阵为Q=[0.01,0;0,0.01],即过程噪声的方差为0.01。取了测量噪声的协方差矩阵为R=[1,0;0,1],即测量噪声的方差为1。根据最小方差估计准则,此时过程噪声方差小于测量噪声的误差,则先验估计值比测量值更可靠。

我们看图:

在速度的分析图中,明显看到速度测量值(灰色)偏离速度真实值(蓝色)的程度大于速度先验估计值(绿色)偏离速度真实值(蓝色)的程度,而经过卡尔曼滤波之后,后验估计值(红色)并没有非常偏离真实值(蓝色)。这就是因为此时卡尔曼滤波更为相信先验估计值。

位置分析同理。

突然想到一个问题:如何确定卡尔曼滤波要迭代多少次呢?

网上说不一定是迭代越多次越准确,由于采用最小方差估计准则,所以我想到了去看误差ek的协方差矩阵的迹,迹越小越好(误差分量的方差之和越小)。然后我又加了几行代码:

# 创建误差的协方差矩阵的迹
    tr_P_posterior = []

# 误差的协方差矩阵的迹,迹越小说明估计越准确
        # print('ek1的方差:',P_posterior[0,0],'ek2的方差',P_posterior[1,1])

        tr_P_posterior.append(P_posterior[0,0]+P_posterior[1,1])

 #误差的协方差的迹
        axs[2].plot(tr_P_posterior,"-",color="blue",label="误差的迹",linewidth="1")

60次迭代的图: 

可以看到,基本上在20次左右,误差的迹就已经收敛至min值了。

于是我把迭代次数调整成20次:

可以看到,大约在十几次的时候,误差的迹就收敛至极限值(约为0.38左右)

那么就是说,刚开始迭代时,卡尔曼滤波器的误差还是挺大的(方差之和大约为1) ,随着迭代的进行,滤波器误差逐步减少至最低点,此后的误差维持在这个点(误差无法完全消除,只存在最小误差),即预测精度达到最优值。

总结一下

1.算法迭代的五个步骤

2.算法的python代码实现

我自己从头开始写的时候最难受的点应该就是因为太久不碰后端,逻辑上会很卡,然后忘记了append函数的作用,搞得我一直在纠结怎么从k-1时刻更新到k时刻,在想是不是要对矩阵做下标的更新什么的,循环迭代这里卡了很久。还有就是对状态变量、先验估计量、后验估计量、协方差矩阵的先验估计量和后验估计量以及它们之间的关系、它们与时刻k、k-1之间的关系不熟。

其实在代码中,我们看一下这五个公式,对于当前时刻k:

step1中的(k-1)时刻的后验估计就是上一次step4估计得到的结果,它们是同一个变量X_posterior;

step2中的 Pk-1就是上一次step5计算得到的结果,它们是同一个变量P_posterior;

step3中的Pk先验,就是本次的step2计算得到的结果,它们是同一个变量P_prior;

step4中的Xk先验,就是本次的step1计算得到的结果,它们是同一个变量X_prior;

其次就是,我们要画图表示出速度、位置的迭代变化,就需要记录下每一次迭代产生的速度值和变量值,然后对它们进行可视化。

最后就是,如果你对先验、后验、时刻搞不清,用英文写清楚变量意思!!不要光贪图简洁!

有关滤波笔记一:卡尔曼滤波(Kalman Filtering)详解的更多相关文章

  1. LC滤波器设计学习笔记(一)滤波电路入门 - 2

    目录前言滤波电路科普主要分类实际情况单位的概念常用评价参数函数型滤波器简单分析滤波电路构成低通滤波器RC低通滤波器RL低通滤波器高通滤波器RC高通滤波器RL高通滤波器部分摘自《LC滤波器设计与制作》,侵权删。前言最近需要学习放大电路和滤波电路,但是由于只在之前做音乐频谱分析仪的时候简单了解过一点点运放,所以也是相当从零开始学习了。滤波电路科普主要分类滤波器:主要是从不同频率的成分中提取出特定频率的信号。有源滤波器:由RC元件与运算放大器组成的滤波器。可滤除某一次或多次谐波,最普通易于采用的无源滤波器结构是将电感与电容串联,可对主要次谐波(3、5、7)构成低阻抗旁路。无源滤波器:无源滤波器,又称

  2. Unity Shader 学习笔记(5)Shader变体、Shader属性定义技巧、自定义材质面板 - 2

    写在之前Shader变体、Shader属性定义技巧、自定义材质面板,这三个知识点任何一个单拿出来都是一套知识体系,不能一概而论,本文章目的在于将学习和实际工作中遇见的问题进行总结,类似于网络笔记之用,方便后续回顾查看,如有以偏概全、不祥不尽之处,还望海涵。1、Shader变体先看一段代码......Properties{ [KeywordEnum(on,off)]USL_USE_COL("IsUseColorMixTex?",int)=0 [Toggle(IS_RED_ON)]_IsRed("IsRed?",int)=0}......//中间省略,后续会有完整代码 #pragmamulti_c

  3. 物联网MQTT协议详解 - 2

    一、什么是MQTT协议MessageQueuingTelemetryTransport:消息队列遥测传输协议。是一种基于客户端-服务端的发布/订阅模式。与HTTP一样,基于TCP/IP协议之上的通讯协议,提供有序、无损、双向连接,由IBM(蓝色巨人)发布。原理:(1)MQTT协议身份和消息格式有三种身份:发布者(Publish)、代理(Broker)(服务器)、订阅者(Subscribe)。其中,消息的发布者和订阅者都是客户端,消息代理是服务器,消息发布者可以同时是订阅者。MQTT传输的消息分为:主题(Topic)和负载(payload)两部分Topic,可以理解为消息的类型,订阅者订阅(Su

  4. Tcl脚本入门笔记详解(一) - 2

    TCL脚本语言简介•TCL(ToolCommandLanguage)是一种解释执行的脚本语言(ScriptingLanguage),它提供了通用的编程能力:支持变量、过程和控制结构;同时TCL还拥有一个功能强大的固有的核心命令集。TCL经常被用于快速原型开发,脚本编程,GUI和测试等方面。•实际上包含了两个部分:一个语言和一个库。首先,Tcl是一种简单的脚本语言,主要使用于发布命令给一些互交程序如文本编辑器、调试器和shell。由于TCL的解释器是用C\C++语言的过程库实现的,因此在某种意义上我们又可以把TCL看作C库,这个库中有丰富的用于扩展TCL命令的C\C++过程和函数,所以,Tcl是

  5. 【详解】Docker安装Elasticsearch7.16.1集群 - 2

    开门见山|拉取镜像dockerpullelasticsearch:7.16.1|配置存放的目录#存放配置文件的文件夹mkdir-p/opt/docker/elasticsearch/node-1/config#存放数据的文件夹mkdir-p/opt/docker/elasticsearch/node-1/data#存放运行日志的文件夹mkdir-p/opt/docker/elasticsearch/node-1/log#存放IK分词插件的文件夹mkdir-p/opt/docker/elasticsearch/node-1/plugins若你使用了moba,直接右键新建即可如上图所示依次类推创建

  6. 【Elasticsearch基础】Elasticsearch索引、文档以及映射操作详解 - 2

    文章目录概念索引相关操作创建索引更新副本查看索引删除索引索引的打开与关闭收缩索引索引别名查询索引别名文档相关操作新建文档查询文档更新文档删除文档映射相关操作查询文档映射创建静态映射创建索引并添加映射概念es中有三个概念要清楚,分别为索引、映射和文档(不用死记硬背,大概有个印象就可以)索引可理解为MySQL数据库;映射可理解为MySQL的表结构;文档可理解为MySQL表中的每行数据静态映射和动态映射上面已经介绍了,映射可理解为MySQL的表结构,在MySQL中,向表中插入数据是需要先创建表结构的;但在es中不必这样,可以直接插入文档,es可以根据插入的文档(数据),动态的创建映射(表结构),这就

  7. 计算机网络笔记:TCP三次握手和四次挥手过程 - 2

    TCP是面向连接的协议,连接的建立和释放是每一次面向连接的通信中必不可少的过程。TCP连接的管理就是使连接的建立和释放都能正常地进行。三次握手TCP连接的建立—三次握手建立TCP连接①若主机A中运行了一个客户进程,当它需要主机B的服务时,就发起TCP连接请求,并在所发送的分段中用SYN=1表示连接请求,并产生一个随机发送序号x,如果连接成功,A将以x作为其发送序号的初始值:seq=x。主机B收到A的连接请求报文,就完成了第一次握手。客户端发送SYN=1表示连接请求客户端发送一个随机发送序号x,如果连接成功,A将以x作为其发送序号的初始值:seq=x②主机B如果同意建立连接,则向主机A发送确认报

  8. 最强Http缓存策略之强缓存和协商缓存的详解与应用实例 - 2

    HTTP缓存是指浏览器或者代理服务器将已经请求过的资源保存到本地,以便下次请求时能够直接从缓存中获取资源,从而减少网络请求次数,提高网页的加载速度和用户体验。缓存分为强缓存和协商缓存两种模式。一.强缓存强缓存是指浏览器直接从本地缓存中获取资源,而不需要向web服务器发出网络请求。这是因为浏览器在第一次请求资源时,服务器会在响应头中添加相关缓存的响应头,以表明该资源的缓存策略。常见的强缓存响应头如下所述:Cache-ControlCache-Control响应头是用于控制强制缓存和协商缓存的缓存策略。该响应头中的指令如下:max-age:指定该资源在本地缓存的最长有效时间,以秒为单位。例如:Ca

  9. IDEA 2022 创建 Spring Boot 项目详解 - 2

    如何用IDEA2022创建并初始化一个SpringBoot项目?目录如何用IDEA2022创建并初始化一个SpringBoot项目?0. 环境说明1.  创建SpringBoot项目 2.编写初始化代码0. 环境说明IDEA2022.3.1JDK1.8SpringBoot1.  创建SpringBoot项目        打开IDEA,选择NewProject创建项目。        填写项目名称、项目构建方式、jdk版本,按需要修改项目文件路径等信息。        选择springboot版本以及需要的包,此处只选择了springweb。        此处需特别注意,若你使用的是jdk1

  10. 华为数通笔记VXLAN&BGP EVPN - 2

    VXLAN简介定义RFC定义了VLAN扩展方案VXLAN(VirtualeXtensibleLocalAreaNetwork,虚拟扩展局域网)。VXLAN采用MACinUDP(UserDatagramProtocol)封装方式,是NVO3(NetworkVirtualizationoverLayer3)中的一种网络虚拟化技术。目的随着网络技术的发展,云计算凭借其在系统利用率高、人力/管理成本低、灵活性/可扩展性强等方面表现出的优势,已经成为目前企业IT建设的新趋势。而服务器虚拟化作为云计算的核心技术之一,得到了越来越多的应用。服务器虚拟化技术的广泛部署,极大地增加了数据中心的计算密度;同时,为

随机推荐