卡尔曼滤波(Kalman Filter)原理浅析-数学理论推导-1

目录

    • 前言
    • 数学理论推导
    • 1. 递归算法
    • 2. 数学基础
    • 结语
    • 参考

前言

最近项目需求涉及到目标跟踪部分,准备从 DeepSORT 多目标跟踪算法入手。DeepSORT 中涉及的内容有点多,以前也就对其进行了简单的了解,但是真正去做发现总是存在这样或者那样的困惑,头疼,果然欠下的总该还的😂

一个个来吧,这个系列文章主要分享博主在学习 DeepSORT 中的 Kalman Filter 的相关知识,主要从两方面分享,一方面是数学理论推导,另一方面是比较通俗易懂的图例分析。

这篇文章主要分享从数学理论推导的方式去理解卡尔曼滤波器,包括递归算法和数学基础

博主为初学者,欢迎交流讨论,若有问题欢迎各位看官批评指正!!!😄

数学理论推导

视频链接:【卡尔曼滤波器】_Kalman_Filter_全网最详细数学推导

:博主也就把 DR_CAN 老师讲解的内容复述了一遍,强烈建议大家观看原视频!!!

1. 递归算法

我们将系统的学习卡尔曼滤波器,由浅入深

首先我们来初尝下递归算法的味道,谈到卡尔曼滤波器(Kalman Filter),这个 Filter 这个词在这里面并不是很能准确地体现出来它的特性。

Kalman Filter 如果用一句话来总结,它是一种 optimal recursive data processing algorithm

它是一种算法,是最优化的、递归的数字处理的算法。它更像是一种观测器,而不是一般意义上的滤波器。卡尔曼滤波器的应用非常广泛尤其是在导航当中,它的广泛应用是因为我们生活的世界中存在着大量的不确定性。当我们去描述一个系统的时候,这个不确定性主要体现在三个方面:

1. 不存在完美的数学模型

2. 系统的扰动往往是不可控的,也是很难建模的

3. 测量的传感器本身存在误差


下面我们来看一个例子,找不同的人用一个尺子去测量一个硬币的直径

在这里插入图片描述

这里面我们用 Z k Z_k Zk 来表示测量的结果,下标 k k k 就代表了第 k k k 次的测量结果,因为每个人测量的值都不一样,而且这个尺子本身也有误差,所以说每次测量的结果也会不尽相同。比如说前三次测量的结果分别为 Z 1 = 50.1 m m Z_1 = 50.1mm Z1=50.1mm Z 2 = 50.4 m m Z_2 = 50.4mm Z2=50.4mm Z 3 = 50.2 m m Z_3 = 50.2mm Z3=50.2mm,这个时候如果我让你去估计一下这个真实结果, 很自然的, 你就会想到去取一个平均值就可以了。

如果用数学来表达,这边我们用 X ^ k \hat{X}_k X^k 来表示第 k k k 次的估计值,其计算如下所示:
X ^ k = 1 k ( Z 1 + Z 2 + ⋅ ⋅ ⋅ + Z k ) = 1 k ( Z 1 + Z 2 + ⋅ ⋅ ⋅ + Z k − 1 ) + 1 k Z k = 1 k k − 1 k − 1 ( Z 1 + Z 2 + ⋅ ⋅ ⋅ + Z k − 1 ) + 1 k Z k = k − 1 k X ^ k − 1 + 1 k Z k = X ^ k − 1 − 1 k X ^ k − 1 + 1 k Z k \begin{aligned} \hat{X}_k &= \frac{1}{k} (Z_1 + Z_2 + \cdot\cdot\cdot + Z_k) \\ &= \frac{1}{k} (Z_1 + Z_2 + \cdot\cdot\cdot + Z_{k-1}) + \frac{1}{k}Z_k \\ &= \frac{1}{k} \frac{k-1}{\color{blue}k-1} \color{blue} (Z_1 + Z_2 + \cdot\cdot\cdot + Z_{k-1})\color{black} + \frac{1}{k}Z_k \\ &= \frac{k-1}{k} \color{blue}\hat{X}_{k-1} \color{black} + \frac{1}{k}Z_k \\ &= \hat{X}_{k-1}-\frac{1}{k}\hat{X}_{k-1}+\frac{1}{k}Z_k \end{aligned} X^k=k1(Z1+Z2++Zk)=k1(Z1+Z2++Zk1)+k1Zk=k1k1k1(Z1+Z2++Zk1)+k1Zk=kk1X^k1+k1Zk=X^k1k1X^k1+k1Zk
重新整理下可得到第 k k k 次数的估计值 X ^ k \hat{X}_k X^k 等于下式:
X ^ k = X ^ k − 1 + 1 k ( Z k − X ^ k − 1 ) \hat{X}_k = \hat{X}_{k-1} + \frac{1}{k}(Z_k - \hat{X}_{k-1}) X^k=X^k1+k1(ZkX^k1)
去分析下这个公式,你可以看到随着 k k k 的增加,即随着次数的增加, 1 k → 0 \frac{1}{k} \to 0 k10 X ^ k → X ^ k − 1 \hat{X}_k \to \hat{X}_{k-1} X^kX^k1,也就是说随着 k k k 的增加,测量的结果就不再重要了。这个也非常好理解,当你有了大量的数据之后,测量了很多次以后,你就对估计的结果很有信心了,所有以后的测量值就变得不那么重要了。而相反当 k k k 值比较小的时候,这个 1 k \frac{1}{k} k1 就会比较大,那这个时候测量的结果就可以起到很大的作用,尤其是它和估计值差距比较大的时候。

下面我们把这个公式提炼出来,我们令 1 k \frac{1}{k} k1 等于 K k \color{blue}K_k Kk,则:
X ^ k = X ^ k − 1 + K k ( Z k − X ^ k − 1 ) \hat{X}_k = \hat{X}_{k-1} + \color{blue}{K}_k \color{black}(Z_k - \hat{X}_{k-1}) X^k=X^k1+Kk(ZkX^k1)
它所代表的意思就是:
当前的估计值 = 上一次的估计值 + 系数 × ( 当前测量值 − 上一次的估计值 ) 当前的估计值 \color{green}= \color{black}上一次的估计值 \color{green} + \color{blue}系数 \color{green} \times \color{black}(当前测量值 \color{green}- \color{black}上一次的估计值) 当前的估计值=上一次的估计值+系数×(当前测量值上一次的估计值)
在卡尔曼滤波器里面,这个 K k \color{blue}K_k Kk 就代表了卡尔曼增益(Kalman Gain),或者叫做卡尔曼因数。通过这个公式可以看出来,新的估计值 X ^ k \hat{X}_k X^k 与上一次的估计值 X ^ k − 1 \hat{X}_{k-1} X^k1 有关,然后上一次的又会与上上次的有关,这就是一种递归的思想。而且这也是卡尔曼滤波器的一个优势,它不需要去追溯很久以前的数据,只需要上一次的就可以了。

而关于这个 Kalman Gain K k \color{blue}{K_k} Kk,在这里先做一个最简单的讨论,让大家有一个初步的认识,首先先引入两个参数,一个是估计误差,也就是估计值和真实值的差距,我们用 e E S T e_{EST} eEST 表示,上面的 e e e 代表 Error 误差,下标的 E S T EST EST 就代表 Estimate 估计,还有一个是测量误差,也就是测量值和真实值的差距,用 e M E A e_{MEA} eMEA 来表示,下标的 M E A MEA MEA 就代表 Measurement 测量

Kalman Gain 就等于下式:
K k = e E S T k − 1 e E S T k − 1 + e M E A k \color{blue}{K_k} \color{black} = \frac{e_{EST_{k-1}}}{e_{EST_{k-1}}+e_{MEA_{k}}} Kk=eESTk1+eMEAkeESTk1
也就是第 k − 1 k-1 k1 次的估计误差除以第 k − 1 k-1 k1 的估计误差加上第 k k k 次的测量误差,这个公式实际上就是卡尔曼滤波器的核心公式。我们会在后面的分析中详细的讲解

它是怎么被推导出来的呢?🤔

这个时候我们先来分析一下它,在 k k k 时刻当 e E S T k − 1 ≫ e M E A k {e_{EST_{k-1}} \gg e_{MEA_{k}}} eESTk1eMEAk K k → 1 \color{blue}{K_k} \to 1 Kk1 X ^ k = X ^ k − 1 + Z k − X ^ k − 1 = Z k \hat{X}_k = \hat{X}_{k-1} + Z_k - \hat{X}_{k-1} = Z_k X^k=X^k1+ZkX^k1=Zk,这就说明了当 k − 1 k-1 k1 时刻的估计误差远大于第 k k k 次的测量误差的时候,这时第 k k k 次的估计值就很趋近于测量值 Z k Z_k Zk 了。这个非常好理解,因为你估计的误差大,测量得更加准确,所以就需要去更加信任这个测量值。

而同理,在 k k k 时刻当 e E S T k − 1 ≪ e M E A k {e_{EST_{k-1}} \ll e_{MEA_{k}}} eESTk1eMEAk K k → 0 \color{blue}{K_k} \to 0 Kk0 X ^ k = X ^ k − 1 \hat{X}_k = \hat{X}_{k-1} X^k=X^k1,也就是说当你的测量误差很大的时候,我们去选择相信这个估计值,相信上一次的估计值。


有了上面的知识,我们先来看一个例子,用一个非常简单的卡尔曼滤波器的思想去解决一个问题,它可以分为三步:

Step1:计算 Kalman Gain K k = e E S T k − 1 e E S T k − 1 + e M E A k \color{blue}K_k \color{black}= \frac{e_{EST_{k-1}}}{e_{EST_{k-1}}+e_{MEA_{k}}} Kk=eESTk1+eMEAkeESTk1

Step2:计算 X ^ k = X ^ k − 1 + K k ( Z k − X ^ k − 1 ) \hat{X}_k = \hat{X}_{k-1} + \color{blue}{K}_k \color{black}(Z_k - \hat{X}_{k-1}) X^k=X^k1+Kk(ZkX^k1)

Step3:更新 e E S T k = ( 1 − K k ) e E S T k − 1 e_{EST_k} = (1 - \color{blue} K_k\color{black})e_{EST_{k-1}} eESTk=(1Kk)eESTk1

有了这三步以后,我们就可以去解决一些实际的问题了。比如说还是测量一个物体它的实际长度,比如说是 50mm,假设我们第一次的估计值 X ^ 0 = 40 m m \hat{X}_0 = 40mm X^0=40mm,第一次的估计误差 e E S T 0 = 5 m m e_{EST_0} = 5mm eEST0=5mm,第一次的测量值 Z 1 = 51 m m Z_1 = 51mm Z1=51mm,测量误差是 e M E A k = 3 m m e_{MEA_k} = 3mm eMEAk=3mm,由于使用同样的一个测量工具,所以它第 k k k 次的测量误差也是 3mm

现在可以做一个表格,然后把有的内容填在上面,我们就可以进行递归的计算了

k k k Z k Z_k Zk e M E A k e_{MEA_k} eMEAk X ^ k \hat{X}_k X^k K k K_k Kk e E S T k e_{EST_k} eESTk
0 40 \color{sienna}40 40 5 \color{red}5 5
1 51 \color{orange}51 51 3 \color{red}3 3

k = 1 k = 1 k=1 时:
K k = 5 5 + 3 = 0.625 X ^ k = 40 + 0.625 ( 51 − 40 ) = 46.875 e E S T = ( 1 − 0.625 ) 5 = 1.875 \begin{aligned} K_k &= \color{red}\frac{5}{5+3} = \color{green}0.625 \\ \hat{X}_k &= \color{sienna}40 + \color{green}0.625(\color{orange}51 -\color{sienna}40) \color{black}= 46.875 \\ e_{EST} &= (1 - \color{green}0.625)\color{red}5 = 1.875 \end{aligned} KkX^keEST=5+35=0.625=40+0.625(5140)=46.875=(10.625)5=1.875
更新下表格:

k k k Z k Z_k Zk e M E A k e_{MEA_k} eMEAk X ^ k \hat{X}_k X^k K k K_k Kk e E S T k e_{EST_k} eESTk
0 40 \color{sienna}40 40 5 \color{red}5 5
1 51 \color{orange}51 51 3 \color{red}3 3 46.875 46.875 46.875 0.625 \color{green}0.625 0.625 1.875 1.875 1.875

k = 2 k = 2 k=2 时:

K k = 1.875 1.875 + 3 = 0.3846 X ^ k = 46.875 + 0.3846 ( 48 − 46.875 ) = 47.308 e E S T = ( 1 − 0.3846 ) 1.875 = 1.154 \begin{aligned} K_k &= \frac{1.875}{1.875+3} = 0.3846 \\ \hat{X}_k &= 46.875 + 0.3846(48-46.875) = 47.308 \\ e_{EST} &= (1 - 0.3846)1.875 = 1.154 \end{aligned} KkX^keEST=1.875+31.875=0.3846=46.875+0.3846(4846.875)=47.308=(10.3846)1.875=1.154
填入到表格中:

k k k Z k Z_k Zk e M E A k e_{MEA_k} eMEAk X ^ k \hat{X}_k X^k K k K_k Kk e E S T k e_{EST_k} eESTk
0 40 \color{sienna}40 40 5 \color{red}5 5
1 51 \color{orange}51 51 3 \color{red}3 3 46.875 46.875 46.875 0.625 \color{green}0.625 0.625 1.875 1.875 1.875
2 48 48 48 3 3 3 47.308 47.308 47.308 0.3846 0.3846 0.3846 1.154 1.154 1.154

后面的递归工作我们交给 Python 来做,代码如下:

import matplotlib.pyplot as plte_MEAk = 3
e_ESTk = 5
Xk = [40] + [0]*16
Zk = [51, 48, 47, 52, 51, 48, 49, 53, 48, 49, 52, 53, 51, 52, 49, 50]Kk = 0
for idx in range(len(Zk)):Kk = e_ESTk / (e_ESTk + e_MEAk)Xk[idx + 1] = Xk[idx] + Kk * (Zk[idx] - Xk[idx])e_ESTk = (1 - Kk) * e_ESTkprint(f"第 {idx:2} 次: 卡尔曼增益 = {Kk:.5f}, 估计值 = {Xk[idx + 1]:.5f}, 估计误差 = {e_ESTk:.5f}")k = [i for i in range(17)]
plt.plot(k, Xk, marker='o', color='red', label='Estimates')
plt.plot(k[:-1], Zk, marker = 'o', color='blue', label='Measurements')
plt.xlabel('Times')
plt.ylabel('Value')
plt.title('Kalman Filter Example')
plt.legend(loc='right')plt.savefig("kf.png", bbox_inches='tight', dpi=300)
plt.show()

输出如下:

0 次: 卡尔曼增益 = 0.62500, 估计值 = 46.87500, 估计误差 = 1.875001 次: 卡尔曼增益 = 0.38462, 估计值 = 47.30769, 估计误差 = 1.153852 次: 卡尔曼增益 = 0.27778, 估计值 = 47.22222, 估计误差 = 0.833333 次: 卡尔曼增益 = 0.21739, 估计值 = 48.26087, 估计误差 = 0.652174 次: 卡尔曼增益 = 0.17857, 估计值 = 48.75000, 估计误差 = 0.535715 次: 卡尔曼增益 = 0.15152, 估计值 = 48.63636, 估计误差 = 0.454556 次: 卡尔曼增益 = 0.13158, 估计值 = 48.68421, 估计误差 = 0.394747 次: 卡尔曼增益 = 0.11628, 估计值 = 49.18605, 估计误差 = 0.348848 次: 卡尔曼增益 = 0.10417, 估计值 = 49.06250, 估计误差 = 0.312509 次: 卡尔曼增益 = 0.09434, 估计值 = 49.05660, 估计误差 = 0.2830210 次: 卡尔曼增益 = 0.08621, 估计值 = 49.31034, 估计误差 = 0.2586211 次: 卡尔曼增益 = 0.07937, 估计值 = 49.60317, 估计误差 = 0.2381012 次: 卡尔曼增益 = 0.07353, 估计值 = 49.70588, 估计误差 = 0.2205913 次: 卡尔曼增益 = 0.06849, 估计值 = 49.86301, 估计误差 = 0.2054814 次: 卡尔曼增益 = 0.06410, 估计值 = 49.80769, 估计误差 = 0.1923115 次: 卡尔曼增益 = 0.06024, 估计值 = 49.81928, 估计误差 = 0.18072

可视化图如下:

在这里插入图片描述

蓝色的是我们的测量结果,红色的是我们一个的估计结果。可以看到我们最开始是从 40 开始的,就是我们初始的估计值,它大概用了五步就到了 49 这个位置,我们之前有说实际值是 50,然后再经过几步的迭代之后,它就越来越靠近我们的实际值。这就是卡尔曼滤波器的一种递归的思想,随着我们不断的去更新,不断的有新的数据进来,它会越来越靠近真实的值。

可以看得出来,这是一个非常简单的利用卡尔曼滤波器的递归思想来做估计的一个例子,我们会在后面的分析中详细展开卡尔曼滤波器的推导过程并给出一些更深层次的应用。

2. 数学基础

这节我们主要讨论四个方面,四个数学基础分别是数据融合(Data Fusion)、协方差矩阵(Covariance Matrix)、状态空间方程(State Space)以及观测器(Observation)

我们先来讨论 Data Fusion 数据融合,还是从一个例子开始,比如说我们用两个称去称一个东西得到了两个结果,分别是 Z 1 = 30 g \color{blue}Z_1=30g Z1=30g Z 2 = 32 g \color{red}Z_2=32g Z2=32g。而且我们知道这两个称都不准,测量都有误差,比如说对于第一个称来说,它的标准差 σ 1 = 2 g \color{blue}\sigma_1=2g σ1=2g 第二个称的标准差 σ 2 = 4 g \color{red}\sigma_2=4g σ2=4g

它们都符合正态分布,也叫高斯分布,所以如果说给出这样的一个条件的话,它们在图中表现出来的就是中型的曲线,比如第一个结果的概率分布如下图:

在这里插入图片描述

中间的位置是 30g,然后标准差是 2,向左向右各一个标准差就覆盖了 68.4% 的可能性,也就是说测出来的结果在 28g~32g 之间的概率是 68.4%

而对于另外一个,因为它的标准差是 4,所以它的图像更矮一些更胖一些,如下图所示:

在这里插入图片描述

关于正态分布和标准差更多细节看观看:【工程数学基础】9_阈值如何选取??在机器视觉中应用正态分布和6-Sigma

我们来把这两个分布绘制在用一张图中,如下图所示:

在这里插入图片描述

这个时候我们有了这两个结果,如果让你去根据这两个结果去估算一下,真实值将会是多少呢?🤔

各位看官可以自己思考一下,如果凭感觉的话,它应该是在这两个结果的中间,而且由于第一个称的误差小(即标准差小),所以它会靠近第一个称称出来的结果。

而如果从数学上去找到一个最优的一个估计值我们又该如何做呢?🤔

这就要用到我们上一节介绍过的 Kalman Gain 的这个思想了,可以令这个估计值等于下式:
Z ^ = Z 1 + K ( Z 2 − Z 1 ) \hat{Z} = \color{blue}Z_1 \color{black} + \color{green}K\color{black}(\color{red}Z_2 \color{black}- \color{blue}Z_1\color{black}) Z^=Z1+K(Z2Z1)
这里面的 K \color{green}K K 就是 Kalman Gain,其值是一个从 0~1 的一个数,可以看出来:

K = 0 K = 0 K=0 时, Z ^ = Z 1 \hat{Z} = \color{blue}Z_1 Z^=Z1

K = 1 K = 1 K=1 时, Z ^ = Z 2 \hat{Z} = \color{red}Z_2 Z^=Z2

下面就是要去求解 K K K 了,我们的目标就是去求解合适的 K K K 使得 Z ^ \hat{Z} Z^ 这个估计值的标准差最小,也就是方差最小。它的方差计算如下:
σ Z ^ 2 = V a r ( Z ^ ) = V a r ( Z 1 + K ( Z 2 − Z 1 ) ) = V a r ( Z 1 − K Z 1 + K Z 2 ) = V a r ( ( 1 − K ) Z 1 + K Z 2 ) = V a r ( ( 1 − K ) Z 1 ) + V a r ( K Z 2 ) = ( 1 − K ) 2 V a r ( Z 1 ) + K 2 V a r ( Z 2 ) = ( 1 − K ) 2 σ 1 2 + K 2 σ 2 2 \begin{aligned} \sigma^2_{\hat{Z}} &= Var(\hat{Z}) \\ &= Var(Z_1 + K(Z_2-Z_1)) \\ &= Var(Z_1 - KZ_1 + KZ_2) \\ &= Var((1-K)Z_1 + KZ_2) \\ &= Var((1-K)Z_1) + Var(KZ_2) \\ &= (1-K)^2Var(Z_1) + K^2Var(Z_2) \\ &= (1-K)^2\sigma_1^2 + K^2\sigma_2^2 \\ \end{aligned} σZ^2=Var(Z^)=Var(Z1+K(Z2Z1))=Var(Z1KZ1+KZ2)=Var((1K)Z1+KZ2)=Var((1K)Z1)+Var(KZ2)=(1K)2Var(Z1)+K2Var(Z2)=(1K)2σ12+K2σ22
其中 Z 1 Z_1 Z1 Z 2 Z_2 Z2 相互独立,我们的目标是要求 σ Z ^ 2 \sigma^2_{\hat{Z}} σZ^2 最小,因此我们可以让它对 K K K 求导,求取其对应的极值,如下:
d σ Z ^ 2 d 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 \begin{aligned} \frac{d\sigma^2_{\hat{Z}}}{dK} &= 0 \\ -2(1-K)\sigma_1^2 + 2K\sigma_2^2 &=0 \\ -\sigma_1^2 + K\sigma_1^2 + K\sigma_2^2 &= 0 \\ K(\sigma_1^2+\sigma_2^2) &= \sigma_1^2 \\ \end{aligned} dKdσZ^22(1K)σ12+2Kσ22σ12+Kσ12+Kσ22K(σ12+σ22)=0=0=0=σ12
因此最终的 K K K 计算如下:
K = σ 1 2 σ 1 2 + σ 2 2 K = \frac{\sigma_1^2}{\sigma_1^2+\sigma_2^2} K=σ12+σ22σ12
我们把 σ 1 = 2 \sigma_1 = 2 σ1=2 σ 2 = 4 \sigma_2 = 4 σ2=4 代入有
K = σ 1 2 σ 1 2 + σ 2 2 = 2 2 2 2 + 4 2 = 4 4 + 16 = 0.2 K = \frac{\sigma_1^2}{\sigma_1^2+\sigma_2^2} = \frac{2^2}{2^2+4^2}=\frac{4}{4+16}=0.2 K=σ12+σ22σ12=22+4222=4+164=0.2
我们把 K = 0.2 K = 0.2 K=0.2 代入有
Z ^ = Z 1 + K ( Z 2 − Z 1 ) = 30 + 0.2 ( 32 − 30 ) = 30.4 \begin{aligned} \hat{Z} &= \color{blue}Z_1 \color{black} + \color{green}K\color{black}(\color{red}Z_2 \color{black}- \color{blue}Z_1\color{black}) \\ &= \color{blue}30 \color{black} + 0.2 (\color{red}32 \color{black}- \color{blue}30\color{black}) \\ &= 30.4 \end{aligned} Z^=Z1+K(Z2Z1)=30+0.2(3230)=30.4
也就是说根据这两个称的特性和测量出来的结果,我们做出了预测,这个预测是 30.4g,并且通过数学证明了它是最优解。这时候也可以去计算它的方差:
σ Z ^ 2 = ( 1 − K ) 2 σ 1 2 + K 2 σ 2 2 = ( 1 − 0.2 ) 2 × 2 2 + 0. 2 2 × 4 2 = 3.2 \begin{aligned} \sigma^2_{\hat{Z}} &= (1-K)^2\sigma_1^2 + K^2\sigma_2^2 \\ &= (1-0.2)^2\times2^2+0.2^2\times4^2 \\ &= 3.2 \end{aligned} σZ^2=(1K)2σ12+K2σ22=(10.2)2×22+0.22×42=3.2
其标准差 σ Z ^ = 3.2 = 1.79 \sigma_{\hat{Z}} = \sqrt{3.2} = 1.79 σZ^=3.2 =1.79

所以说我们的这个最优的估计值 σ Z ^ 2 = 30.4 g \sigma^2_{\hat{Z}} = 30.4g σZ^2=30.4g 标准差 σ Z ^ = 1.79 \sigma_{\hat{Z}} = 1.79 σZ^=1.79

如果将它绘制在之前的图中,它应该是一个更高更瘦的图像。如下图所示:

在这里插入图片描述

这个过程就叫做数据融合,大家要掌握这个思想,记住这个理念,过一会我们会再回过头来看它。


第二个要聊的就是协方差矩阵 Covariance Matrix,它把方差和协方差在一个矩阵当中表现出来,也就是体现了变量间的一个联动关系。还是从一个例子入手

球员身高体重年龄
瓦尔迪1797433
奥巴梅扬1878031
萨拉赫1757128
Mean180.37530.7

上表是 2020 年英超中射手榜前三球员的一些基本的数据,包括身高、体重和年龄。我们将它们当作三个变量分别为 x x x y y y z z z,然后我们求下它们的平均值:
x ˉ = 1 3 ( 179 + 187 + 175 ) = 180.3 y ˉ = 1 3 ( 74 + 80 + 71 ) = 75 z ˉ = 1 3 ( 33 + 31 + 28 ) = 30.7 \begin{aligned} \bar x &= \frac{1}{3}(179+187+175) = 180.3 \\ \bar y &= \frac{1}{3}(74+80+71) = 75 \\ \bar z &= \frac{1}{3}(33+31+28) = 30.7 \end{aligned} xˉyˉzˉ=31(179+187+175)=180.3=31(74+80+71)=75=31(33+31+28)=30.7
分别是身高 180.3,体重 75,年龄 30.7,然后我们还可以计算下它们对应的方差:
σ x 2 = 1 3 ( ( 179 − 180.3 ) 2 + ( 187 − 180.3 ) 2 + ( 175 − 180.3 ) 2 ) = 24.89 σ y 2 = 1 3 ( ( 74 − 75 ) 2 + ( 80 − 75 ) 2 + ( 71 − 75 ) 2 ) = 14 σ z 2 = 1 3 ( ( 33 − 30.7 ) 2 + ( 31 − 30.7 ) 2 + ( 28 − 30.7 ) 2 ) = 4.22 \begin{aligned} \sigma_x^2 &= \frac{1}{3}((179-180.3)^2+(187-180.3)^2+(175-180.3)^2) = 24.89 \\ \sigma_y^2 &= \frac{1}{3}((74-75)^2+(80-75)^2+(71-75)^2) = 14 \\ \sigma_z^2 &= \frac{1}{3}((33-30.7)^2+(31-30.7)^2+(28-30.7)^2) = 4.22 \end{aligned} σx2σy2σz2=31((179180.3)2+(187180.3)2+(175180.3)2)=24.89=31((7475)2+(8075)2+(7175)2)=14=31((3330.7)2+(3130.7)2+(2830.7)2)=4.22
协方差计算如下:
σ x σ y = 1 3 ( ( 179 − 180.3 ) ( 74 − 75 ) + ( 187 − 180.3 ) ( 80 − 75 ) + ( 175 − 180.3 ) ( 71 − 75 ) ) = 18.7 = σ y σ x σ x σ z = 1 3 ( ( 179 − 180.3 ) ( 33 − 30.7 ) + ( 187 − 180.3 ) ( 31 − 30.7 ) + ( 175 − 180.3 ) ( 28 − 30.7 ) ) = 4.4 = σ z σ x σ y σ z = 1 3 ( ( 74 − 75 ) ( 33 − 30.7 ) + ( 80 − 75 ) ( 31 − 30.7 ) + ( 71 − 75 ) ( 28 − 30.7 ) ) = 3.3 = σ z σ y \begin{aligned} \sigma_x \sigma_y &= \frac{1}{3}((179-180.3)(74-75)+(187-180.3)(80-75)+(175-180.3)(71-75)) = 18.7 = \sigma_y \sigma_x\\ \sigma_x \sigma_z &= \frac{1}{3}((179-180.3)(33-30.7)+(187-180.3)(31-30.7)+(175-180.3)(28-30.7)) = 4.4 = \sigma_z \sigma_x \\ \sigma_y \sigma_z &= \frac{1}{3}((74-75)(33-30.7)+(80-75)(31-30.7)+(71-75)(28-30.7)) = 3.3 = \sigma_z \sigma_y \end{aligned} σxσyσxσzσyσz=31((179180.3)(7475)+(187180.3)(8075)+(175180.3)(7175))=18.7=σyσx=31((179180.3)(3330.7)+(187180.3)(3130.7)+(175180.3)(2830.7))=4.4=σzσx=31((7475)(3330.7)+(8075)(3130.7)+(7175)(2830.7))=3.3=σzσy
可以看到如果最后计算出来的值大于 0,则说明两个变量的变化方向是一样的,这时候我们称这两个变量正相关。而如果计算出来的值小于 0,则说明两个变量的变化方向是相反的,这时候我们称这两个变量负相关

然后说到协方差矩阵 P P P,它的形式如下所示:
P = [ σ x 2 σ x σ y σ x σ z σ y σ x σ y 2 σ y σ z σ z σ x σ z σ y σ z 2 ] = [ 24.89 18.7 4.4 18.7 14 3.3 4.4 3.3 4.22 ] P=\begin{bmatrix} \sigma_x^2 & \sigma_x\sigma_y & \sigma_x\sigma_z \\ \sigma_y\sigma_x & \sigma_y^2 & \sigma_y\sigma_z \\ \sigma_z\sigma_x & \sigma_z\sigma_y & \sigma_z^2 \\ \end{bmatrix} =\begin{bmatrix} 24.89 & 18.7 & 4.4 \\ 18.7 & 14 & 3.3 \\ 4.4 & 3.3 & 4.22 \\ \end{bmatrix} P= σx2σyσxσzσxσxσyσy2σzσyσxσzσyσzσz2 = 24.8918.74.418.7143.34.43.34.22
然后我们就可以分析各个数据之间的关系,因为这里只取了三组数据,所以不是非常准确,我们就希望把它扩展一下,在扩展之前我们先来看看如何通过矩阵的运算来计算这个协方差。

这个对于以后要编程的话非常有帮助,我们还是看这个例子,首先我们要先求一个过渡矩阵 a a a

a = [ x 1 y 1 z 1 x 2 y 2 z 2 x 3 y 3 z 3 ] − 1 3 [ 1 1 1 1 1 1 1 1 1 ] [ x 1 y 1 z 1 x 2 y 2 z 2 x 3 y 3 z 3 ] a=\begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ \end{bmatrix}- \frac{1}{3} \begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \\ \end{bmatrix} \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ \end{bmatrix} a= x1x2x3y1y2y3z1z2z3 31 111111111 x1x2x3y1y2y3z1z2z3
协方差矩阵 P P P 的计算则为 P = 1 3 a T a P = \frac{1}{3}a^Ta P=31aTa,通过这样的方法,我们就可以算出协方差矩阵了。

我们现在来分析一下这个协方差矩阵说明了什么,下表中给出了 15 位球员的一个基本信息,它涵盖了五大联赛在 2020 年射手榜的前三位:

球员身高体重年龄
瓦尔迪1797433
奥巴梅扬1878031
萨拉赫1757128
梅西1707233
本泽马1858132
杰拉德1777528
因莫比莱1787130
C罗1878335
卢卡库1909427
姆巴佩1787321
本耶德尔1706829
邓贝莱1837323
莱万多夫斯基1847831
维尔纳1807524
桑乔1807630
Mean180.276.2728.3

下表则是它们的协方差矩阵:

身高体重年龄
身高32.6929.751.4
体重29.7538.063.98
年龄1.43.97819.4

大家可以仔细观察下,先看协方差矩阵对角线的这几个数,可以看到身高体重的变化还是蛮大的,它们的方差还是很大的。这也就说明了,如果你想成为射手并不会很局限于身高和体重,年龄的跨度也比较大,说明一个射手和年龄的关系也不是特别的大。

接着我们再来看下它们的协方差的关系,首先是身高和体重,可以看到它们的协方差等于 29.75,也就是说明体重和身高是非常的正相关的,随着身高的增加,体重也增加,这是非常非常合理的。然后我们再看年龄和体重和身高,它们的相关性就非常小,一个是 1.4 一个是 3.98,这就说明了对于这些运动员来说,他们的身高、体重和年龄就没有太大的一个关系了,就说明他们确实是保持的还是很好的,并没有因为年龄大了就变胖


最后一个话题就是状态空间的表达式和观测器,其实关于状态空间的表达来说,它是一个完整的体系。现代控制理论就是以状态空间方程为基础的,感兴趣的可参考视频:现代控制理论

在这里我们只是给出一些最基本的一些概念,依然从一个例子入手,一个典型的弹簧振动阻尼系统,如下图所示:

在这里插入图片描述

它的质量是 m m m,向右施加的力是 F F F,向右的方向是 x x x,弹簧的胡可系数是 k k k,阻尼系数是 B B B,它的动态方程的表达式如下:
m x ¨ + B x ˙ + k x = F m\ddot{x}+B\dot{x}+kx=F mx¨+Bx˙+kx=F
我们可以把 F F F 定义成 u u u 也就是系统的输入,这个时候如果要把它化成状态空间的一种表达形式,首先就要确定两个状态变量,我们可以令 x 1 = x x_1 = x x1=x x 2 = x ˙ x_2 = \dot{x} x2=x˙,这样的话有
x ˙ 1 = x 2 x ˙ 2 = x ¨ = 1 m u − B m x ˙ − k m x = 1 m u − B m x 2 − k m x 1 \begin{aligned} \dot{x}_1 &= x_2 \\ \dot{x}_2 &= \ddot{x} \\ &= \frac{1}{m}u-\frac{B}{m}\dot{x}-\frac{k}{m}x \\ &= \frac{1}{m}u-\frac{B}{m}x_2-\frac{k}{m}x_1 \end{aligned} x˙1x˙2=x2=x¨=m1umBx˙mkx=m1umBx2mkx1
这样就用两个一阶微分方程把这个形式表达出来了,如果我们这里面还有一个测量量,我们可以把它定义成 z 1 z_1 z1 z 1 z_1 z1 测的是它的位置 z 1 = x = x 1 z_1 = x = x_1 z1=x=x1 z 2 z_2 z2 测的是它的速度 z 2 = x ˙ = x 2 z_2 = \dot{x} = x_2 z2=x˙=x2,这样的话就是一个测量的一个方程。

然后我们可以把上面的四个式子用矩阵的方式把它很紧凑的表达出来,如下所示:
[ x ˙ 1 x ˙ 2 ] = [ 0 1 − k m − B m ] [ x 1 x 2 ] + [ 0 1 m ] u [ z 1 z 2 ] = [ 1 0 0 1 ] [ x 1 x 2 ] \begin{aligned} \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \end{bmatrix} &= \begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{B}{m} \\ \end{bmatrix} \begin{bmatrix} {x}_1 \\ {x}_2 \\ \end{bmatrix} + \begin{bmatrix} 0 \\ \frac{1}{m} \\ \end{bmatrix} u \\ \\ \begin{bmatrix} z_1 \\ z_2 \\ \end{bmatrix} &= \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ \end{bmatrix} \begin{bmatrix} {x}_1 \\ {x}_2 \\ \end{bmatrix} \end{aligned} [x˙1x˙2][z1z2]=[0mk1mB][x1x2]+[0m1]u=[1001][x1x2]
这就是状态空间的表达形式,归纳出来如下所示:
X ˙ ( t ) = A X ( t ) + B U ( t ) Z ( t ) = H X ( t ) \begin{aligned} \dot{X}(t) &= AX(t)+BU(t) \\ Z(t) &= HX(t) \end{aligned} X˙(t)Z(t)=AX(t)+BU(t)=HX(t)
上式是一种连续的表达形式,我们看到 X ˙ ( t ) \dot{X}(t) X˙(t) 就是 X X X 对时间 t t t 的导数,它体现了随时间的变化。如果我们把它写成离散的一种形式,每相隔一个时刻来看这个变化的话,它就可以写成下面的形式:
X k = A X k − 1 + B U k Z k = H X k \begin{aligned} X_k &= AX_{k-1}+BU_{k} \\ Z_k &= HX_k \end{aligned} XkZk=AXk1+BUk=HXk
这个下标 k k k k − 1 k-1 k1 k + 1 k+1 k+1,每一个单位 1 都代表了一个时间单位叫做 sample time,也就是采样时间。这样的话它也代表了一种变化趋势。

当然从连续型到离散型的表达不是直接照抄上面的矩阵结果,是需要根据采样时间来计算的,但这部分内容不是这里所讨论的重点,这里面你只需要理解它的基本概念就可以了。它是体现了一种变化,是从上一步到这一步的一种变化。

好,最后如果大家还记得之前所说的,世界中充满了不确定性。这个时候我们就可以加一些不确定性在里边了。比如说我们的模型也许不准确,加一个 W k − 1 \color{red}W_{k-1} Wk1 它叫做过程噪音(Process Noise),关于测量我们也可以加上 V k \color{red}V_k Vk 它叫做测量噪音(Measurement Noise),是在测量当中产生的不确定性,那我们整个表达式就变成了下面的形式:
X k = A X k − 1 + B U k + W k − 1 Z k = H X k + V k \begin{aligned} X_k &= AX_{k-1}+BU_{k}\color{red}+W_{k-1} \\ Z_k &= HX_k\color{red}+V_{k} \end{aligned} XkZk=AXk1+BUk+Wk1=HXk+Vk
现在我们就有了两个不确定性了,也就是说模型不准确,然后测出来的值也不准确,在这两个都不确定的情况下,如何去估计一个精确的 x ^ k \hat{x}_k x^k 呢?🤔

这就是卡尔曼滤波器需要解决的问题!

大家回想一下我们在开头讲解的数据融合的例子,我们现在遇到的情况其实和上面的那个情况也是差不多的,只不过我们现在有一个不太准的计算结果和一个不太准的测量结果,我们要根据这两个结果来估计一个相对准确的值,来找到一个误差比它们两个本身都要小的一个结果。

这就是我们后面要重点去分析的内容,到此为止,关于 Kalman Filter 的基本理念就讲完了,概念性的直观理解也到此结束了,从之后开始就是比较枯燥的数学推导了,请大家做好准备😄。

结语

本篇博客首先从递归算法出发来理解卡尔曼滤波器,卡尔曼滤波器是一种最优化的、递归的数字处理算法,它更像是一种观测器,而不是一般意义上的滤波器。接着我们讲解了一个实际的例子,利用卡尔曼滤波器的递归思想来估计硬币的直径,随着我们不断地去更新,我们估计的值会越来越接近实际值,这是递归思想的体现。

随后我们讲解了推导卡尔曼滤波器需要的数学基础,包括数据融合、协方差矩阵、状态空间方程和观测器四个部分。我们从两个不准确的称称同一个物体如何获得一个比较准确的结果这个例子出发讲解了数据融合的概念,我们从运动员的身高体重年龄之间的相关性这个例子出发讲解了协方差矩阵的概念,最后我们通过弹簧振动阻尼系统讲解了状态空间表达式的概念。

在下一篇文章中我们将会详细推导卡尔曼增益和误差协方差矩阵,敬请期待😄

参考

  • 现代控制理论
  • 常用的MarkDown颜色 笔记
  • 【卡尔曼滤波器】_Kalman_Filter_全网最详细数学推导【博主强烈推荐!!!👍👍👍】
  • 【工程数学基础】9_阈值如何选取??在机器视觉中应用正态分布和6-Sigma

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/80968.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Java idea编译器工程out目录修改

问题 多个工程在一个文件夹下,有时会变为所有的工程只用一个out文件夹,这时运行会出错。 解决方案是将各自的工程out放在自己的文件夹下面 以上面这个工程为例 如果project structure里面有则直接按照下面的指标设置,如果没有则添加到里面再…

ABB 1TGE120010R1300 控制主板模块

ABB 1TGE120010R1300 控制主板模块是一种用于控制和监测电力设备的模块,具有以下功能: 控制和监测电力设备:该模块可以通过与电力设备连接来控制和监测设备的性能和状态,例如启停设备、调节电压和功率等。 通信功能:该…

docker 和k8s 入门

docker 和k8s 入门 本文是云原生的学习记录,可以参考以下文档 k8s https://www.yuque.com/leifengyang/oncloud 相关视频教程可参考如下 https://www.bilibili.com/video/BV13Q4y1C7hS?p2&vd_source0882f549dac54045384d4a921596e234 相对于公有云&#x…

【力扣周赛】第 362 场周赛(⭐差分匹配状态压缩DP矩阵快速幂优化DPKMP)

文章目录 竞赛链接Q1:2848. 与车相交的点解法1——排序后枚举解法2——差分数组⭐差分数组相关题目列表📕1094. 拼车1109. 航班预订统计2381. 字母移位 II2406. 将区间分为最少组数解法1——排序贪心优先队列解法2——差分数组 2772. 使数组中的所有元素…

TCP IP网络编程(五) 基于TCP的服务器端、客户端 (补充)

文章目录 回声客户端的完美实现回声客户端出现的问题回声客户端问题解决方法 TCP原理TCP套接字中的I/O缓冲TCP内部工作原理1:与对方套接字的连接TCP内部工作原理2:与对方主机的数据交换TCP内部工作原理3:断开与套接字的连接 回声客户端的完美…

双向链表的实现(增删查改)——最好理解的链表

双向链表的实现 一,双向链表的特点二,双向链表的结构三,双向链表的内容实现3.1创建node节点3.2初始化3.3打印3.4插入3.4.1尾插3.4.2头插3.4.3在pos位置上插入 3.5删除3.5.1尾删3.5.2头删3.5.3删除pos位置上的数据 四,调试技巧&…

day21算法

常见的七种查找算法: ​ 数据结构是数据存储的方式,算法是数据计算的方式。所以在开发中,算法和数据结构息息相关。今天的讲义中会涉及部分数据结构的专业名词,如果各位铁粉有疑惑,可以先看一下哥们后面录制的数据结构…

【C++学习】继承

目录 一、继承的概念及定义 1、继承的概念 2、继承的定义 2.1 定义格式 2.2 继承关系和访问限定符 2.3 继承基类成员访问方式的变化 二、基类和派生类对象赋值转换 三、继承中的作用域 四、派生类的默认成员函数 五、继承与友元 六、继承与静态成员 七、复杂的菱形…

SpringMvc第五战-【SpringMvcJSR303和拦截器】

前言: 小编阐述了springmvc 中的文件下载,以及jrebel的使用和文件下载以及多文件下载! 在本次小编将会介绍JSR303的概念,应用场景和在具体实例的使用;和拦截器的应用 一.JSR303的介绍 1.什么是JSR303? JSR是Java S…

Aztec的隐私抽象:在尊重EVM合约开发习惯的情况下实现智能合约隐私

1. 引言 Aztec的架构,不同于当前“通过EVM兼容执行环境”所实现的区块链水平扩容趋势。Aztec内部笑称其构建的为首个非zkEVM协议。 Aztec专注于实现: 成为理解和需要智能合约隐私的开发者的终极解决方案。 Aztec为开发者提供构建隐私优先app所需的网…

【微信小程序开发】宠物预约医疗项目实战-环境配置与Vant UI集成

第一章 宠物预约医疗项目实战-环境配置与Vant UI集成 文章目录 前言一、Vant UI是什么?二、使用步骤2.1 安装 node.js2.2 通过 npm 安装vant2.3 修改 app.json2.4 修改 project.config.json2.5 构建 npm 包2.6 使用组件全局引入和局部引入全局引入局部引入 前言 Va…

基于Java+SpringBoot+Vue+uniapp点餐小程序(亮点:协同过滤算法、会员系统,购物车结算、在线聊天)

校园点餐小程序 一、前言二、我的优势2.1 自己的网站2.2 自己的小程序(小蔡coding)2.3 有保障的售后2.4 福利 三、开发环境与技术3.1 MySQL数据库3.2 Vue前端技术3.3 Spring Boot框架3.4 微信小程序 四、功能设计4.1 系统功能结构设计4.2 主要功能描述 五…

SpringBoot整合Flowable

1. 配置 &#xff08;1&#xff09; 引入maven依赖 <dependency><groupId>org.flowable</groupId><artifactId>flowable-spring-boot-starter</artifactId><version>6.7.2</version></dependency><!-- MySQL连接 -->&l…

MySQL--MySQL索引事务

事务的概念 事务指逻辑上的一组操作&#xff0c;组成这组操作的各个单元&#xff0c;要么全部成功&#xff0c;要么全部失败。 在不同的环境中&#xff0c;都可以有事务。对应在数据库中&#xff0c;就是数据库事务。 使用 &#xff08;1&#xff09;开启事务&#xff1a;start…

什么是接口自动化?为什么要做?和怎么做接口自动化?

服务端接口测试介绍 什么是服务端&#xff1f; 一般所说的服务端是指为用户在 APP 或 PC 使用的互联网功能提供数据服务的背后的一切。以天猫精灵智能音箱系列的产品链路为例&#xff0c;服务端便是网关&#xff08;包括网关在内&#xff09;之后的链路。 什么是接口&#xf…

【自然语言处理】【大模型】RWKV:基于RNN的LLM

相关博客 【自然语言处理】【大模型】RWKV&#xff1a;基于RNN的LLM 【自然语言处理】【大模型】CodeGen&#xff1a;一个用于多轮程序合成的代码大语言模型 【自然语言处理】【大模型】CodeGeeX&#xff1a;用于代码生成的多语言预训练模型 【自然语言处理】【大模型】LaMDA&a…

深入网络底层,了解Linux系统收发网络数据包的过程、原理、流程,附图文说明

深入网络底层&#xff0c;了解Linux系统收发网络数据包的过程、原理、流程&#xff0c;附图文说明。 Linux 服务器收到网络数据包&#xff0c;需要经过哪些处理&#xff0c;一步步将数据传给应用进程的呢&#xff1f;应用进程发送数据包时&#xff0c;Linux 又是如何操作将数据…

android studio platform使用体验分享(as无法跳转c/c++等native源码的福音,强烈推荐)

hi&#xff0c;粉丝朋友们&#xff1a; 大家好&#xff01;这些天粉丝朋友们分享了一下Android Studio for Platform 这个最新的google开发的阅读aosp源码的工具&#xff0c;特别适合做原生系统开发。具体官方介绍如下地址&#xff1a; 参考链接&#xff1a;https://developer.…

react的状态管理简单钩子方法

1.recoil useProvider文件: import { atom, useRecoilState } from recoil;const initState atom({key: initState,default: {state: [],}, })// 将业务逻辑拆分到一个单独文件中&#xff0c;方便进行状态管理 export interface StateProps {id: number;text: string;isFini…

异地远程访问本地SQL Server数据库【无公网IP内网穿透】

文章目录 1. 前言2. SeaFile云盘设置2.1 Owncould的安装环境设置2.2 SeaFile下载安装2.3 SeaFile的配置 3. cpolar内网穿透3.1 Cpolar下载安装3.2 Cpolar的注册3.3 Cpolar云端设置3.4 Cpolar本地设置 4. 公网访问测试5. 结语 1. 前言 现在我们身边的只能设备越来越多&#xff…