草庐IT

快速上手的Python版二维卡尔曼滤波解析

佐咖 2023-04-03 原文

卡尔曼滤波是最好的线性滤波,但是需要推导的公式教多,也很细,这里推荐一个B站博主视频讲解的关于卡尔曼滤波,讲的很好,很细,适合小白学习,链接地址为:添加链接描述。如果完全没接触过卡尔曼滤波的,建议从第一集开始学习。
下面是我跟着这位博主学习后,再加上其他大神写的代码,融入我自己的理解,对代码进行修改后的版本,每一个部分都有详细的注释,更加的通俗易懂,希望能帮助到需要快速上手卡尔曼滤波的学习者。

卡尔曼的本质就下面五个公式:

参数说明,见下:

在卡尔曼滤波中需要调的参数主要有两个:
(1)Q矩阵:过程的协方差矩阵。注意:过程的协方差越小,说明预估的更准确,更可信。

(2)R矩阵:测量的协方差矩阵。注意:测量的协方差矩阵越小,说明测量的结果更准确,更可信。

怎么理解呢,假设你的测量设备是几万元或几千元的,测量值的可靠度较高(测量值的权值占比就需要较大),那么你的R矩阵就定义的小一些,比如定义为[[0.1,0],[0,0.1],这时候Q矩阵就可以定义的大一些,比如定义为[[1,0],[0,1]]。

如果你是使用的测量设备是某宝9.9元包邮的,那么你就应该相信估计值(估计值的权重占比就需要较大),这时候你的R矩阵就定义的大一些,比如定义为[[1,0],[0,1],这时候Q矩阵就可以定义的小一些,比如定义为[[0.11,0],[0,0.11]]。

想要获得较好的滤波预测效果,得不断的调试这两个矩阵的大小。

下面是详细的代码部分:

import numpy as np
import matplotlib.pyplot as plt

##下面是通过高斯正态分布产生误差
def gaussian_distribution_generator(var):
    return np.random.normal(loc=0.0, scale=var, size=None)  ##np.random.normal是正态分布

# 参数loc(float):正态分布的均值,对应着这个分布的中心。loc=0说明这一个以Y轴为对称轴的正态分布,
# 参数scale(float):正态分布的标准差,对应分布的宽度,scale越大,正态分布的曲线越矮胖,scale越小,曲线越高瘦。
# 参数size(int 或者整数元组):输出的值赋在shape里,默认为None。

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

# 过程噪声协方差矩阵Q,p(w)~N(0,Q),噪声来自真实世界中的不确定性      N(0,Q) 表示期望是0,协方差矩阵是Q
Q = np.array([[0.1, 0],              # Q其实就是过程的协方差
              [0, 0.1]])             # 过程的协方差越小,说明预估的越准确

# 观测噪声协方差矩阵R,p(v)~N(0,R)     也是测量的协方差矩阵
R = np.array([[1, 0],                # R其实就是测量的协方差矩阵
              [0, 1]])               #测量的协方差越小,说明测量的结果更准确

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

# 控制输入矩阵B
B = None

# 初始位置与速度
X0 = np.array([[0],                  # 二维的,定义初始的位置和速度,初始的位置为0 ,速度为1
               [1]])

# 状态估计协方差矩阵P初始化(其实就是初始化最优解的协方差矩阵)
P = np.array([[1, 0],
              [0, 1]])

if __name__ == "__main__":
    # ---------------------------初始化-------------------------
    X_true = np.array(X0)           # 真实状态初始化
    X_posterior = np.array(X0)      # X_posterior表示上一时刻的最优估计值
    P_posterior = np.array(P)       # P_posterior是继续更新最优解的协方差矩阵

    speed_true = []                 # 实际的速度值
    position_true = []              # 实际的位置值

    speed_measure = []             # 测量到的速度值
    position_measure = []          # 测量到的位置值

    speed_prior_est = []           # 速度的先验误差 估计值
    position_prior_est = []        # 位置的先验误差 估计值

    speed_posterior_est = []       # 根据估计值及当前时刻的观测值融合到一体得到的最优估计速度值存入到列表中
    position_posterior_est = []    # 根据估计值及当前时刻的观测值融合到一体得到的最优估计值位置值存入到列表中

    for i in range(30):
        # -----------------------生成真实值----------------------
        # 生成过程噪声
        # 通过Q矩阵产生误差
        w = np.array(
            [[gaussian_distribution_generator(Q[0, 0])],         # gaussian_distribution_generator(Q[0, 0])的值为0.01103097596,(这值是随机的,每次运行都不一样)
             [gaussian_distribution_generator(Q[1, 1])]])        # gaussian_distribution_generator(Q[1, 1])的值为-0.1242726240,(这值是随机的,每次运行都不一样)
        X_true = np.dot(A, X_true) + w                           # 得到当前时刻实际的速度值和位置值,其中A是状态转移矩阵上一时刻的状态转移到当前时刻
        speed_true.append(X_true[1, 0])                          # 将第二行第一列位置的数值,1.01103098添加到列表speed_true里面
        position_true.append(X_true[0, 0])                       # 将第一行第一列位置的数值,1.01103098添加到列表position_true里面

        # -----------------------生成观测值----------------------
        # 生成观测噪声
        ##通过R矩阵产生误差
        v = np.array(
            [[gaussian_distribution_generator(R[0, 0])],         # gaussian_distribution_generator(R[0, 0])的值为-0.62251186549,(这值是随机的,每次运行都不一样)
             [gaussian_distribution_generator(R[1, 1])]])        # gaussian_distribution_generator(R[1, 1])的值为2.52779100481,(这值是随机的,每次运行都不一样)
        Z_measure = np.dot(H, X_true) + v                        # 生成观测值,H为单位阵E    # A是状态观测矩阵    # Z_measure表示测量得到的值
        position_measure.append(Z_measure[0, 0])
        speed_measure.append(Z_measure[1, 0])

        ################################################################下面开始进行预测和更新,来回不断的迭代#######################################################################

        # ----------------------进行先验估计---------------------
        X_prior = np.dot(A,X_posterior)                          # X_prior表示根据上一时刻的最优估计值得到当前的估计值  X_posterior表示上一时刻的最优估计值
        position_prior_est.append(X_prior[0, 0])                 # 将根据上一时刻计算得到的最优估计位置值添加到列表position_prior_est中
        speed_prior_est.append(X_prior[1, 0])                    # 将根据上一时刻计算得到的最优估计速度值添加到列表speed_prior_est.append中

        # 计算状态估计协方差矩阵P
        P_prior_1 = np.dot(A,P_posterior)                        # P_posterior是上一时刻最优估计的协方差矩阵   # P_prior_1就为公式中的(A.Pk-1)
        P_prior = np.dot(P_prior_1, A.T) + Q                     # P_prior是得出当前的先验估计协方差矩阵    #Q是过程协方差

        # ----------------------计算卡尔曼增益,用numpy一步一步计算Prior and posterior
        k1 = np.dot(P_prior, H.T)                                # P_prior是得出当前的先验估计协方差矩阵
        k2 = np.dot(np.dot(H, P_prior), H.T) + R                 # R是测量的协方差矩阵
        K = np.dot(k1, np.linalg.inv(k2))                        # np.linalg.inv():矩阵求逆   # K就是当前时刻的卡尔曼增益

        # ---------------------后验估计------------
        X_posterior_1 = Z_measure - np.dot(H, X_prior)           # X_prior表示根据上一时刻的最优估计值得到当前的估计值   # Z_measure表示测量得到的值
        X_posterior = X_prior + np.dot(K, X_posterior_1)         # K就是当前时刻的卡尔曼增益    # X_posterior是根据估计值及当前时刻的观测值融合到一体得到的最优估计值
        position_posterior_est.append(X_posterior[0, 0])         # 根据估计值及当前时刻的观测值融合到一体得到的最优估计位置值存入到列表中
        speed_posterior_est.append(X_posterior[1, 0])            # 根据估计值及当前时刻的观测值融合到一体得到的最优估计速度值存入到列表中

        # 更新状态估计协方差矩阵P     (其实就是继续更新最优解的协方差矩阵)
        P_posterior_1 = np.eye(2) - np.dot(K, H)                 # np.eye(2)返回一个二维数组,对角线上为1,其他地方为0
        P_posterior = np.dot(P_posterior_1, P_prior)             # P_posterior是继续更新最优解的协方差矩阵  ##P_prior是得出当前的先验估计协方差矩阵

    # 可视化显示
    if True:
        plt.rcParams['font.sans-serif'] = ['SimHei']  # 坐标图像中显示中文
        plt.rcParams['axes.unicode_minus'] = False

        fig, axs = plt.subplots(1, 2)
        axs[0].plot(speed_true, "-", label="speed_true(实际值)", linewidth=1)  # Plot some data on the axes.
        axs[0].plot(speed_measure, "-", label="speed_measure(测量值)", linewidth=1)  # Plot some data on the axes.
        axs[0].plot(speed_prior_est, "-", label="speed_prior_est(估计值)", linewidth=1)  # Plot some data on the axes.
        axs[0].plot(speed_posterior_est, "-", label="speed_posterior_est(融合测量值和估计值)",
                    linewidth=1)  # Plot some data on the axes.
        axs[0].set_title("speed")
        axs[0].set_xlabel('k')  # Add an x-label to the axes.
        axs[0].legend()  # Add a legend.

        axs[1].plot(position_true, "-", label="position_true(实际值)", linewidth=1)  # Plot some data on the axes.
        axs[1].plot(position_measure, "-", label="position_measure(测量值)", linewidth=1)  # Plot some data on the axes.
        axs[1].plot(position_prior_est, "-", label="position_prior_est(估计值)",
                    linewidth=1)  # Plot some data on the axes.
        axs[1].plot(position_posterior_est, "-", label="position_posterior_est(融合测量值和估计值)",
                    linewidth=1)  # Plot some data on the axes.
        axs[1].set_title("position")
        axs[1].set_xlabel('k')  # Add an x-label to the axes.
        axs[1].legend()  # Add a legend.

        plt.show()

最后的运行结果见下:

以上就是二维卡尔曼滤波的相关解析和代码,仔细看,仔细品,你看得懂的。

有关快速上手的Python版二维卡尔曼滤波解析的更多相关文章

  1. Ruby 解析字符串 - 2

    我有一个字符串input="maybe(thisis|thatwas)some((nice|ugly)(day|night)|(strange(weather|time)))"Ruby中解析该字符串的最佳方法是什么?我的意思是脚本应该能够像这样构建句子:maybethisissomeuglynightmaybethatwassomenicenightmaybethiswassomestrangetime等等,你明白了......我应该一个字符一个字符地读取字符串并构建一个带有堆栈的状态机来存储括号值以供以后计算,还是有更好的方法?也许为此目的准备了一个开箱即用的库?

  2. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  3. ruby - 解析 RDFa、微数据等的最佳方式是什么,使用统一的模式/词汇(例如 schema.org)存储和显示信息 - 2

    我主要使用Ruby来执行此操作,但到目前为止我的攻击计划如下:使用gemsrdf、rdf-rdfa和rdf-microdata或mida来解析给定任何URI的数据。我认为最好映射到像schema.org这样的统一模式,例如使用这个yaml文件,它试图描述数据词汇表和opengraph到schema.org之间的转换:#SchemaXtoschema.orgconversion#data-vocabularyDV:name:namestreet-address:streetAddressregion:addressRegionlocality:addressLocalityphoto:i

  4. ruby - 用逗号、双引号和编码解析 csv - 2

    我正在使用ruby​​1.9解析以下带有MacRoman字符的csv文件#encoding:ISO-8859-1#csv_parse.csvName,main-dialogue"Marceu","Giveittohimóhe,hiswife."我做了以下解析。require'csv'input_string=File.read("../csv_parse.rb").force_encoding("ISO-8859-1").encode("UTF-8")#=>"Name,main-dialogue\r\n\"Marceu\",\"Giveittohim\x97he,hiswife.\"\

  5. ruby-on-rails - 我更新了 ruby​​ gems,现在到处都收到解析树错误和弃用警告! - 2

    简而言之错误:NOTE:Gem::SourceIndex#add_specisdeprecated,useSpecification.add_spec.Itwillberemovedonorafter2011-11-01.Gem::SourceIndex#add_speccalledfrom/opt/local/lib/ruby/site_ruby/1.8/rubygems/source_index.rb:91./opt/local/lib/ruby/gems/1.8/gems/rails-2.3.8/lib/rails/gem_dependency.rb:275:in`==':und

  6. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  7. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  8. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  9. python - 如何读取 MIDI 文件、更改其乐器并将其写回? - 2

    我想解析一个已经存在的.mid文件,改变它的乐器,例如从“acousticgrandpiano”到“violin”,然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容,该乐器通过program_change或patch_change指令进行了更改,但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作,但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成,每个音轨是十六个channel中任何一个上的

  10. 「Python|Selenium|场景案例」如何定位iframe中的元素? - 2

    本文主要介绍在使用Selenium进行自动化测试或者任务时,对于使用了iframe的页面,如何定位iframe中的元素文章目录场景描述解决方案具体代码场景描述当我们在使用Selenium进行自动化测试的时候,可能会遇到一些界面或者窗体是使用HTML的iframe标签进行承载的。对于iframe中的标签,如果直接查找是无法找到的,会抛出没有找到元素的异常。比如近在咫尺的例子就是,CSDN的登录窗体就是使用的iframe,大家可以尝试通过F12开发者模式查看到的tag_name,class_name,id或者xpath来定位中的页面元素,会抛出NoSuchElementException异常。解决

随机推荐