卡尔曼滤波(Kalman Filtering)——(3)数据融合 状态空间方程

article/2024/12/23 17:17:01

一、数据融合

假设举例

  假设测量一物体的质量,现在有两个测量设备但是都存在误差且误差服从正态分布,分别用两个设备进行测量第一设备测量的测量值为 z 1 z_1 z1 ,标准差为 σ 1 \sigma_1 σ1;第二设备测量的测量值为 z 2 z_2 z2 ,标准差为 σ 2 \sigma_2 σ2 。现在有两个测量值那么应该更愿意去相信那个设备测出的值呢?或者把两次数据进行处理导出一个更接近真实值的值呢?
  如果已知: z 1 = 30 g ; σ 1 = 2 g z_1=30g;\sigma_1=2g z1=30gσ1=2g
        z 2 = 32 g ; σ 2 = 4 g z_2=32g;\sigma_2=4g z2=32gσ2=4g
  现在让估算真实值,该如何计算呢?
分析:显而易见最终所估算的真实值也是服从正态分布的,且其结果在 z 1 z_1 z1 z 2 z_2 z2之间,由于标准差 σ 1 比 σ 2 \sigma_1比\sigma_2 σ1σ2小则更靠近 z 1 z_1 z1

公式推导过程

  通过上一博客可知:当前的估计值 = 上一次的估计值 + 系数 (当前测量值 - 上一次的估计值)Kalman Filtering通过调整观测值和估计值的权重,使修正后的值更趋向于相信观测值还是相信估计值的一个过程。当我们使用同一设备进行测量时估计值通常是进行数学建模计算得出,由于有的误差无法建模所以得出的值只能是估计值,如第一篇博客所举的例子中当前估计值与上一次估计值分别是第 K K K次取平均值和第 K − 1 K-1 K1次取平均值,而当前测量值为第 K K K次测量结果。从这我们也可以看出数据融合的思想,数学建模值与测量值之间融合得到最接近真实值的结果。但是本博客举得例子是两个设备测量结果没有进行建模,应该如何使用Kalman Filtering?其实同理只需要将一个值作为数学建模的结果即估计值,另一个作为测量结果即可。所以Kalman Filtering的数据融合是两个值之间的融合,不论值是建模得来还是传感器测得都行。
  现在要得出一个最优估计值 z ^ \hat{z} z^,由以上面推导的公式可知:(将 z 1 z_1 z1当作估计值; z 2 z_2 z2当作测量值)
z ^ = z 1 + K ( z 2 − z 1 ) (1) \hat{z}=z_{1}+K\left(z_{2}-z_{1}\right)\tag1 z^=z1+K(z2z1)(1)
式中 : K ⟶ : K \longrightarrow :K 卡尔曼增益/因数, 且 K ∈ [ 0 , 1 ] K \in[0,1] K[0,1]
   K = 0 K=0 K=0 ⇒ \Rightarrow z ^ = z 1 \hat{z}=z_{1} z^=z1
   K = 1 K=1 K=1 ⇒ \Rightarrow z ^ = z 2 \hat{z}=z_{2} z^=z2
现在目标是求 K K K 使得 σ z ^ \sigma_{\hat{z}} σz^ 最小 ⇒ \Rightarrow K K K 使得方差 Var ⁡ ( z ^ ) \operatorname{Var}(\hat{z}) Var(z^) 最小 ⇒ \Rightarrow K K K使误差最小
则:
σ z ^ 2 = V ar ⁡ ( z 1 + K ( z 2 − z 1 ) ) = V ar ⁡ ( z 1 − K z 1 + K z 2 ) = V ar ⁡ ( ( 1 − K ) z 1 + K z 2 ) ⟶ 两 个 设 备 测 量 结 果 与 误 差 , 相 互 独 立 , 互 不 影 响 = Var ⁡ ( ( 1 − K ) z 1 ) + Var ⁡ ( K z 2 ) ⟶ K 为 常 数 , 从 方 差 中 提 出 得 带 平 方 = ( 1 − K ) 2 Var ⁡ ( z 1 ) + K 2 Var ⁡ ( z 2 ) = ( 1 − K ) 2 σ 1 2 + K 2 σ 2 2 \begin{aligned} \sigma_{\hat{z}}^{2} &=\boldsymbol{V} \operatorname{ar}\left(z_{1}+K\left(z_{2}-z_{1}\right)\right) \\ &=\boldsymbol{V} \operatorname{ar}\left(z_{1}-K z_{1}+K z_{2}\right) \\ &=\boldsymbol{V} \operatorname{ar}\left((1-K) z_{1}+K z_{2}\right) \longrightarrow两个设备测量结果与误差,相互独立,互不影响\\ &=\operatorname{Var}\left((1-K) z_{1}\right)+\operatorname{Var}\left(K z_{2}\right) \longrightarrow K为常数,从方差中提出得带平方\\ &=(1-K)^{2} \operatorname{Var}\left(z_{1}\right)+K^{2} \operatorname{Var}\left(z_{2}\right) \\ &=(1-K)^{2} \sigma_{1}^{2}+K^{2} \sigma_{2}^{2} \end{aligned} σz^2=Var(z1+K(z2z1))=Var(z1Kz1+Kz2)=Var((1K)z1+Kz2)=Var((1K)z1)+Var(Kz2)K=(1K)2Var(z1)+K2Var(z2)=(1K)2σ12+K2σ22

求方差最小值时,则方差 σ ı ^ 2 \sigma_{\hat{\imath}}^{2} σı^2 K K K 的导数为零,即:
∂ σ z ^ 2 ∂ K = 0 ⇒ − 2 ( 1 − K ) σ 1 2 + 2 K σ 2 2 = 0 ⇒ − σ 1 2 + K σ 1 2 + K σ 2 2 = 0 ⇒ K ( σ 1 2 + σ 2 2 ) = σ 1 2 ⇒ K = σ 1 2 σ 1 2 + σ 2 2 \begin{aligned} \frac{\partial \sigma_{\hat{z}}^{2}}{\partial K}=0 & \Rightarrow-2(1-K) \sigma_{1}^{2}+2 K \sigma_{2}^{2}=0 \\ & \Rightarrow-\sigma_{1}^2+K \sigma_{1}^{2}+K \sigma_{2}^{2}=0 \\ & \Rightarrow K\left(\sigma_{1}^{2}+\sigma_{2}^{2}\right)=\sigma_{1}^{2} \\ & \Rightarrow K=\frac{\sigma_{1}^{2}}{\sigma_{1}^{2}+\sigma_{2}^{2}} \end{aligned} Kσz^2=02(1K)σ12+2Kσ22=0σ12+Kσ12+Kσ22=0K(σ12+σ22)=σ12K=σ12+σ22σ12
可以看出K的取值与进行数据融合的两个值的误差大小有关:
   当 σ 2 > > σ 1 \sigma_{2}>>\sigma_{1} σ2>>σ1 K = 0 K=0 K=0 ⇒ \Rightarrow z ^ = z 1 \hat{z}=z_{1} z^=z1
   当 σ 2 < < σ 1 \sigma_{2}<<\sigma_{1} σ2<<σ1 K = 1 K=1 K=1 ⇒ \Rightarrow z ^ = z 2 \hat{z}=z_{2} z^=z2
σ 1 = 2 g ; σ 2 = 4 g \sigma_1=2g;\sigma_2=4g σ1=2g;σ2=4g带入数据得: K = 0.2 K=0.2 K=0.2带入(1)式中得 z ^ = 30.4 g \hat{z}=30.4g z^=30.4g
       σ z ^ 2 = 3.2 ⇒ σ z ^ = 1.79 < m i n [ σ 1 2 , σ 2 2 ] \sigma_{\hat{z}}^{2}=3.2 \Rightarrow \sigma_{\hat{z}}=1.79<min{[\sigma_1^2,\sigma_2^2]} σz^2=3.2σz^=1.79<min[σ12,σ22]
上面所进行的过程就是数据融合的过程,由此可见将两个设备所测得的数据进行融合所得到的误差更小更加接近真实值。卡尔曼滤波利用信息融合减少误差从而使估计值更加接近真实值也就是寻找最优估计得过程。

再次理解

  现在可以看出卡尔曼滤波信息融合的关键就是要求卡尔曼增益 K K K,这个就是一个权值使当前估计值趋近于误差小的估计值。上面通过严谨的数学公式推导得出了估计值 z ^ = 30.4 g \hat{z}=30.4g z^=30.4g,下面用小学二年级的知识在来理解理解:
  已知: z 1 = 30 g ; σ 1 = 2 g z_1=30g;\sigma_1=2g z1=30gσ1=2g
      z 2 = 32 g ; σ 2 = 4 g z_2=32g;\sigma_2=4g z2=32gσ2=4g
解:
    加权平均值
z ^ = z 1 σ 2 2 σ 1 2 + σ 2 2 + z 2 σ 1 2 σ 1 2 + σ 2 2 = 30 × 4 2 2 2 + 4 2 + 32 × 2 2 2 2 + 4 2 = 30.4 \begin{aligned} \hat{z}&=z_1\frac{\sigma _2^2}{\sigma _1^2+\sigma _2^2} +z_2\frac{\sigma _1^2}{\sigma _1^2+\sigma _2^2}\\ &=30\times\frac{4^2}{2^2+4^2} + 32\times\frac{2^2}{2^2+4^2} \\ &=30.4 \end{aligned} z^=z1σ12+σ22σ22+z2σ12+σ22σ12=30×22+4242+32×22+4222=30.4
  两个设备测量的结果不相关,相互独立,则两个的加权平均值就为 z ^ \hat{z} z^。其中误差小的测量结果乘以权重大的;反之误差大的测量结果乘以权重小的;很好理解因为我们是比较相信误差小的结果的。
  由于两个测量误差不相关则:(可以理解为并联电阻
1 σ z ^ 2 = 1 σ 1 2 + 1 σ 2 2 = 1 2 2 + 1 4 2 = 0.3125 ⇒ σ z ^ 2 = 3.2 ⇒ σ z ^ = 1.79 \begin{aligned} \frac{1}{\sigma _{\hat{z} }^2} &=\frac{1}{\sigma _{1 }^2}+\frac{1}{\sigma _{2 }^2}\\ &=\frac{1}{2^2}+\frac{1}{4^2}\\ &=0.3125\\ \Rightarrow \sigma _{\hat{z} }^2&=3.2\Rightarrow\sigma _{\hat{z} }=1.79 \end{aligned} σz^21σz^2=σ121+σ221=221+421=0.3125=3.2σz^=1.79
如果误差相关误差相关则为: σ z ^ 2 = σ 1 2 + σ 2 2 \sigma _{\hat{z}}^2 = \sigma _{{1}}^2 + \sigma _{{2} }^2 σz^2=σ12+σ22(串联电阻)
  一般情况下:串联误差为误差相关
           并联误差为误差不相关

二、状态空间方程

  系统会有一系列的状态,这个状态你可以任意的去选择,只要你觉得能够充分表达你所需要的信息。系统状态会发生变化,那么这个变化应当和什么有关呢?

  1. 系统当前的状态;
  2. 系统的输入。

阻尼滑块模型

  一个弹簧滑块系统,滑块的加速度(受到的外力)和系统当前的位置有关(弹簧弹力),也和外力输入有关(人手或电机给与滑块一个外力)。将输入记为u。系统状态的变化,对于连续系统而言,是系统的导数,对于离散系统而言,是系统状态的变化量,一般取为系统下一时刻的状态( x k + 1 \boldsymbol{x}_{k+1} xk+1= A x k + B u k \boldsymbol{A} \boldsymbol{x}_{k}+B u_{k} Axk+Buk),二者包含的信息一致。
在这里插入图片描述

系统的输出和什么有关呢?当然是和系统当前状态和输入有关。

1、连续表达式

该系统动态方程表达式为:
m x ¨ + c x ˙ + k x = F m \ddot{x}+c \dot{x}+k x=F mx¨+cx˙+kx=F
x 1 = x , x 2 = x ˙ , u = F x_{1}=x, x_{2}=\dot{x}, u=F x1=x,x2=x˙,u=F, 有
x ˙ 1 = x 2 x ˙ 2 = x ¨ = 1 m u − c m x ˙ − k m x = 1 m u − c m x 2 − k m x 1 \begin{gathered} \dot{x}_{1}=x_{2} \\ \dot{x}_{2}=\ddot{x}=\frac{1}{m} u-\frac{c}{m} \dot{x}-\frac{k}{m} x \\ =\frac{1}{m} u-\frac{c}{m} x_{2}-\frac{k}{m} x_{1} \end{gathered} x˙1=x2x˙2=x¨=m1umcx˙mkx=m1umcx2mkx1
假设位置测量量为 z 1 = x = x 1 z_{1}=x=x_{1} z1=x=x1, 速度测量量为 z 2 = x ˙ = x 2 z_{2}=\dot{x}=x_{2} z2=x˙=x2,将上式用矩阵形式表达则有:
[ x ˙ 1 x ˙ 2 ] = [ 0 1 − k m − c m ] [ x 1 x 2 ] + [ 0 1 m ] u \begin{aligned} \left[\begin{array}{c} \dot{x}_{1} \\ \dot{x}_{2} \end{array}\right]=\left[\begin{array}{cc} 0 & 1 \\ -\frac{k}{m} & -\frac{c}{m} \end{array}\right]\left[\begin{array}{l} x_{1} \\ x_{2} \end{array}\right]+\left[\begin{array}{c} 0 \\ \frac{1}{m} \end{array}\right] u \end{aligned} [x˙1x˙2]=[0mk1mc][x1x2]+[0m1]u
[ z 1 z 2 ] = [ 1 0 0 1 ] [ x 1 x 2 ] \begin{aligned} \left[\begin{array}{l} z_{1} \\ z_{2} \end{array}\right]=\left[\begin{array}{ll} 1 & 0 \\ 0 & 1 \end{array}\right]\left[\begin{array}{l} x_{1} \\ x_{2} \end{array}\right] \end{aligned} [z1z2]=[1001][x1x2]
上述表达形式即状态空间表达形式,上述表达式可分别归纳为随时间变化表达式:
x ˙ ( t ) = A x ( t ) + B u ( t ) z ( t ) = H x ( t ) \begin{gathered} \dot{\boldsymbol{x}}(t)=A \boldsymbol{x}(t)+B u(t) \\ \boldsymbol{z}(t)=H \boldsymbol{x}(t) \end{gathered} x˙(t)=Ax(t)+Bu(t)z(t)=Hx(t)

2、离散表达形式

x k = A x k − 1 + B u k − 1 z k = H x k \begin{gathered} \boldsymbol{x}_{k}=\boldsymbol{A} \boldsymbol{x}_{k-1}+B u_{k-1} \\ \boldsymbol{z}_{k}=H \boldsymbol{x}_{k} \end{gathered} xk=Axk1+Buk1zk=Hxk
其中,下标 k 、 k − 1 、 k + 1 k、 k-1、k+1 kk1k+1 表示采样时间。即系统当前状态与系统上一时刻状态和上一时刻输入有关。

  由于所有的系统数学模型往往是理想模型,存在许多假设,与实际模型存在差距。比如上述阻尼滑块模型表达式中就忽略了阻力的存在即系统存在噪声且噪声的表达式未知,所以应当加入未知噪声时表达式为:
x k = A x k − 1 + B u k − 1 + w k − 1 z k = H x k + v k \begin{gathered} \boldsymbol{x}_{k}=\boldsymbol{A} \boldsymbol{x}_{k-1}+B u_{k-1}+\boldsymbol{w}_{k-1} \\ \boldsymbol{z}_{k}=H \boldsymbol{x}_{k}+\boldsymbol{v}_{k} \end{gathered} xk=Axk1+Buk1+wk1zk=Hxk+vk
其中: w ⟶ w \longrightarrow w 过程噪音
    v ⟶ v \longrightarrow v 测量噪音
卡尔曼滤波器的目的即从存在噪声的数据中求出误差最小的结果

3、符号含义

符号意义备注
x k x_{k} xk系统状态矩阵(真实值)计算得到
z k z_{k} zk状态阵的观测量实际测量
A A A状态转移矩阵
B B B控制输入矩阵
H H H状态观测矩阵
w k − 1 w_{k-1} wk1过程噪声服从正态分布
v k v_{k} vk测量噪声服从正态分布

参考文献

  DR_CAN
  【官方教程】卡尔曼滤波器教程与MATLAB仿真(全)(中英字幕)


http://chatgpt.dhexx.cn/article/lrxTqnL8.shtml

相关文章

[现代控制理论]2_state-space状态空间方程

[现代控制理论]11_现代控制理论串讲_完结_pdf获取 [现代控制理论]10_可观测性与分离原理_观测器与控制器 [现代控制理论]9_状态观测器设计_龙伯格观测器 [现代控制理论]8.5_线性控制器设计_轨迹跟踪simulink [现代控制理论]8_LQR控制器_simulink [现代控制理论]7_线性控制器设计…

现代控制理论(一) 状态空间方程

文章目录 状态方程和输出方程基本的状态空间方程 线性非线性时变时不变系统的能控能观状态方程的解无输入线性时不变转移矩阵拉普拉斯求转移矩阵 有输入线性时不变 控制什么&#xff1f; 输入u1,u2,u3,…输出y1,y2,y3…的系统。u是控制量&#xff0c;y是响应结果&#xff0c;也…

【离散系统】传递函数和状态空间方程离散化

本文如有错误&#xff0c;恳请指正。 目录 离散系统 采样控制系统 数字控制系统 信号采样 采样定理&#xff08;香农定理&#xff09; 信号保持—零阶保持器 Z变换 Z 变换方法 级数求和法 部分分式法 基本定理 Z反变换 Z反变换方法 长除法 部分分式法&#xff0…

matlab 状态空间的波特图,MATLAB:对于状态空间方程的系统辨识

本文介绍了如何利用MATLAB辨识状态空间方程中的未知参数。 假设我们的被控系统的表达如下: 我们想要通过实验数据辨识出参数K1和K2​,方法如下: 第一步,采集实验数据。 需要的数据包括系统一段时间内的系统输出Y(ts),以及控制量U(ts),这些数据应该是以某个固定的采样频率…

状态空间方程系统建模

以质量弹簧阻尼系统为例&#xff0c;它的动态微分方程之前提到过为&#xff0c; 令此系统的输入等于外力&#xff0c;系统的输出等于位移。 现代控制理论使用状态空间方程的表达方式。 状态空间——一个集合&#xff0c;输入、输出及状态变量&#xff0c;用一系列一阶方程表达…

状态空间方程转换传递函数

对状态空间方程公式(1)进行拉氏变换 对状态空间方程公式(2)进行拉氏变换 公式(5)带入公式(3)&#xff0c;得到输出和输入的关系 最终转换为传递函数表示

现代控制工程笔记(一)控制系统的状态空间描述

文章目录 1. 基本概念2. 系统的状态空间描述状态空间描述框图状态变量选取的非唯一性 3. 由系统微分方程列写状态空间表达式一.微分方程中不包含输入函数的导数项相变量法其他方法&#xff1a; 二.微分方程中包含输入函数的导数项 4. 由传递函数列写状态空间表达式直接实现串联…

鲁迅《狂人日记》全文

看了《觉醒年代》&#xff0c;不觉找到鲁迅的《狂人日记》全文&#xff0c;摘录在这里。 狂人日记 鲁迅 狂人日记序 某君昆仲&#xff0c;今隐其名&#xff0c;皆余昔日在中学时良友&#xff1b;分隔多年&#xff0c;消息渐阙。日前偶闻其一大病&#xff1b;适归故乡&#xff0…

趣味三角——第5章——苍穹和地球的测量

目录 5.1. 三角学在测量天体和地球的应用中发展 5.2. Abraham De Moivre在天体和地球测量中的数学贡献 第5章 天体和地球的测量 5.1. 三角学在测量苍穹和地球的应用中发展 The science of trigonometry was in a sense a precursor of the telescope. It brought farawa…

Linux 设备驱动程序(一)

Linux 内核系列文章 Linux 内核设计与实现 深入理解 Linux 内核 Linux 设备驱动程序&#xff08;一&#xff09; Linux 设备驱动程序&#xff08;二&#xff09; Linux 设备驱动程序&#xff08;三&#xff09; Linux 设备驱动程序&#xff08;四&#xff09; Linux设备驱动开发…

【Linux】进程和计划任务管理

文章目录 一、线程、进程、程序的概念什么是程序&#xff1f;什么是线程&#xff1f;什么是进程&#xff1f;线程与进程的关系线程与进程的区别程序与进程的区别 二、查看进程的方式查看静态的进程统计信息——psps -aux命令ps -elf命令ps查看线程命令 过滤查询——grep查看动态…

从“账房先生”到“中国巨型计算机之父”,慈云桂先后主导了中国四代计算机的研发...

作者 | 年素清责编 | 李雪敬出品 | CSDN&#xff08;ID&#xff1a;CSDNnews&#xff09; 慈云桂是中国计算机界的先驱&#xff0c;他主导了中国第一台电子管计算机的研制&#xff0c;接着是晶体管计算机、集成电路计算机和“银河”亿次巨型机&#xff0c;是公认的“中国巨型计…

《钢琴调律原理及应用》 笔记

【【调律理论篇】】 【第一章 绪论】第一节 钢琴调律的概念 美国人威廉布雷德怀特于 1917 年发表了世界上第一部关于钢琴调律理论与技术的著作&#xff0c;书名为《钢琴调律与相关技术》 福岛琢郎于1950年发表一部名为《钢琴的构造调律修理》的专著 80年代初&#xff0c;在沈…

C#文件操作从入门到精通(1)——INI文件操作

前言: 我们在开发c#的winform程序中,因为有些参数是不断变化的 ,所以经常需要开放一些参数提供给用户设置,通过操作Ini文件来保存我们设置的参数也是c#开发中经常使用的技术,本文就来详细介绍操作ini文件的以下功能: 1、读取ini文件,获取某个节点的某个键的值 2、写入i…

【Android Framework系列】第7章 WMS原理

1 前言 前面【Android Framework系列】第5章 AMS启动流程和【Android Framework系列】第6章 AMS原理之Launcher启动流程我们分析了AMS启动以及Launcher启动的整体流程&#xff0c;那Launcher(Activity启动)后&#xff0c;UI是如何渲染到屏幕并且展示出来的呢&#xff1f;我们这…

知识图谱的建立与查询(以党史人物查询为例)

目录 0 前言 1.确定实体关系属性 2.通过EasyDL标注 3.抽取出实体和关系 4.查询 5.前端页面 0 前言 今天终于答完辩了&#xff0c;真的是舌战群儒&#xff0c;好在有惊无险。贴出唯一一张拍的答辩现场照片。&#xff08;照片里没有我哈哈哈&#xff09; 1.确定实体关系属…

[译]理解PG如何执行一个查询-1

理解PG如何执行一个查询 PG服务器收到客户端发来的查询后&#xff0c;查询的文本交给解析器。解析器扫描查询并检查它的语法。若语法正确&#xff0c;解析器会将查询文本转换成解析树。解析树是一种以正式、明确的形式表示查询含义的数据结构。给定查询&#xff1a; SELECT cus…

【*一篇足以*Java并发编程实践】《Java并发编程实践》学习Note - Part3

目录&#xff1a; 1.避免活跃度危险 1.1 死锁 1.2 避免和诊断死锁 1.3.其他活跃度危险 2.性能和可伸缩性 2.1 内存同步 2.2 阻塞 2.3 减少锁的竞争 3.Lock、ReentrantLock和Synchronized 3.1 可轮询和可定时的锁请求 3.1 可中断的锁获取操作 4.原子变量与非阻塞同步…

阿朱说:咨询的历史(万字深度长文)

&#xff08;1&#xff09;知识成为资产&#xff1a;瓦特蒸汽机 13世纪的英国&#xff0c;首先产生了人类历史上的第一部专利保护法。不过最初是很粗糙的&#xff0c;授予专利的权力完全掌握在国王手中&#xff0c;发放专利特许证&#xff0c;将某种独占经营权授予工匠、商人&a…

SpringBoot + Thymeleaf 练手小项目 --------- 豆瓣网站模拟

目录 一、项目介绍二、资源准备1. 准备数据库表2. 准备image、css、js等静态资源文件3. 项目结构 三、开发步骤1. 新建项目2. pom.xml3. 实体类 model4. Mapper 类5. service 类6. 首页 index.html 开发① MovieController② index.html 7. 电影详情页 movie_info.html 开发① …