数据处理的常用方式
在讨论卡尔曼滤波器的思想之前,先从大多数人熟悉的环节开始讲起。在实际的工作生活中,当我们对同一事物进行多次采样后获得一组数据时,往往会对这组数据进行求均值,以此作为这组数据的具有一定优化意义的估计值。例如当我们对一个物品进行称重时,往往会对这一物品进行多次称重记录数据,并用这组数据的平均值来作为这一物品的最终重量等等。这一类多次测量取平均的方法,实际上包含了数据处理中的优化估计的思维在里面,但这种优化方式并不一定是最优的。
另外,在很多情况下,我们都需要实时的对系统中某个状态进行估计,因此我们会在每个采样周期结束后就会进行一次估计,估计的方式同样是用前面采样结果的平均值来作为当前时刻的估计值。这种估计方式会随着采样数量的增加而逐渐收敛于状态的实际值附近。对于一个连续稳定间隔采样的数据集来说,
a
k
(
k
=
1
,
2
,
…
,
n
)
a_{k}(k=1,2, \ldots, n)
ak(k=1,2,…,n)表示在
k
k
k时刻采样得到的数据,通过求前k时刻采样数据的平均值,来得到k时刻的估计值
a
^
k
\hat{a}_{k}
a^k。此时其代数形式为
a
^
k
=
1
k
(
a
1
+
a
2
+
…
+
a
k
)
\hat{a}_{k}=\frac{1}{\mathrm{k}}\left(a_{1}+a_{2}+\ldots+a_{k}\right)
a^k=k1(a1+a2+…+ak)
从上式可以看出,在对每个时刻进行平均估计的时候,需要记录之前每一个时刻的采样数据,当采样进行到一定程度之后,这种方式会带来负担,也即其空间复杂度和时间复杂度都不够合理。这个时候可以引入迭代的思想来解决这个问题。
a
^
k
=
1
k
(
a
1
+
a
2
+
…
+
a
k
)
a
^
k
=
1
k
k
−
1
k
−
1
(
a
1
+
a
2
+
…
+
a
k
)
a
^
k
=
1
k
k
−
1
k
−
1
a
k
+
1
k
k
−
1
k
−
1
(
a
1
+
a
2
+
…
+
a
k
−
1
)
a
^
k
=
1
k
a
k
+
k
−
1
k
1
k
−
1
(
a
1
+
a
2
+
…
+
a
k
−
1
)
a
^
k
=
1
k
a
k
+
k
−
1
k
a
^
k
−
1
a
^
k
=
1
k
a
k
+
(
1
−
1
k
)
a
^
k
−
1
\hat{a}_{k}=\frac{1}{k}\left(a_{1}+a_{2}+\ldots+a_{k}\right) \\ \hat{a}_{k}=\frac{1}{k} \frac{k-1}{k-1}\left(a_{1}+a_{2}+\ldots+a_{k}\right) \\ \hat{a}_{k}=\frac{1}{k} \frac{k-1}{k-1} a_{k}+\frac{1}{k} \frac{k-1}{k-1}\left(a_{1}+a_{2}+\ldots+a_{k-1}\right) \\ \hat{a}_{k}=\frac{1}{k} a_{k}+\frac{k-1}{k} \frac{1}{k-1}\left(a_{1}+a_{2}+\ldots+a_{k-1}\right) \\ \hat{a}_{k}=\frac{1}{k} a_{k}+\frac{k-1}{k} \hat{a}_{k-1} \\ \hat{a}_{k}=\frac{1}{k} a_{k}+\left(1-\frac{1}{k}\right) \hat{a}_{k-1}
a^k=k1(a1+a2+…+ak)a^k=k1k−1k−1(a1+a2+…+ak)a^k=k1k−1k−1ak+k1k−1k−1(a1+a2+…+ak−1)a^k=k1ak+kk−1k−11(a1+a2+…+ak−1)a^k=k1ak+kk−1a^k−1a^k=k1ak+(1−k1)a^k−1
从结果可以看出,随着
k
k
k的增大,
k
k
k时刻的估计值会更受到
k
−
1
k-1
k−1时刻的估计值的影响,直到
k
k
k趋向于无穷时,估计值将收敛。从这里开始,让我联想起之前做数据处理时用到最多的低通滤波器公式。但不同于现在的是这个权重系数往往是根据滤波效果调试出来的固定值,并没有一个计算的过程。
注:原先使用的最多的低通滤波器的公式是
a
^
k
=
K
∗
a
k
+
(
1
−
K
)
a
^
k
−
1
\hat{a}_{k}=K * a_{k}+(1-K) \hat{a}_{k-1}
a^k=K∗ak+(1−K)a^k−1
简单来书哦就是用当前最新的传感器采样数据和上一时刻滤波结果的加权和得到当前时刻的滤波结果,从上述的公式推到中可以看出,这种滤波方法的本质还是最开始提到的平均估计。
这里涉及到一个问题,我们应该如何去权衡上一时刻估计值和当前时刻采样值之间的关系。换句话说,我们应该更相信上一时刻的估计值,还是当前时刻的采样值。这个问题直接决定了加权系数
K
K
K的大小,也决定了估计后的值的趋势。
这个问题本质为数据融合的问题,即当有对同一事物的多个不同数据时,我们如何以最优的方式来把他们都利用起来,而不是不由分说的轻信其中某几组数据。这里的最优也需要有严格的数学证明来保证。
数据融合
下面用一个例子来讨论数据融合。
假设对一个物品进行测重,有两个不一样的称,两个称都不完全准确,存在着偏差。在数学上,一般用方差来描述这种偏差。因此对于两个称的称重结果,我们做出这样的数学表达:
{
z
1
=
5.2
g
,
δ
1
=
0.2
z
2
=
4.9
g
,
δ
2
=
0.4
\left\{\begin{array}{l} z_{1}=5.2 g, \delta_{1}=0.2 \\ z_{2}=4.9 g, \delta_{2}=0.4 \end{array}\right.
{z1=5.2g,δ1=0.2z2=4.9g,δ2=0.4
其中,
z
1
z_{1}
z1和
z
2
z_{2}
z2分别表示两个称的测量结果,而
δ
1
\delta_{1}
δ1和
δ
2
\delta_{2}
δ2表示两个称的方差,代表的物理意义即称本身存在的测量误差。这种误差导致的不准确性可以用高斯分布来描述,根据概率论的相关知识,可以得到一下关系
{
z
1
∼
N
(
5.2
,
0.2
)
z
2
∼
N
(
4.9
,
0.4
)
\left\{\begin{array}{l} z_{1} \sim N(5.2,0.2) \\ z_{2} \sim N(4.9,0.4) \end{array}\right.
{z1∼N(5.2,0.2)z2∼N(4.9,0.4)
第一个称称出来的结果是5.2g,那这个物品实际的重量就满足一个概率分布,即以5.2为数学期望,0.2为方差的高斯分布。具体的概率分布的意义就不在此赘述,自行补充。
可以比较简单粗暴的理解为,第一个称得到的结果是5.2g,那物品的实际重量大概率会在5~5.4之间,当然也有可能是其他的重量,只不过概率相对较小而已。
根据之前的描述,我们现在要在这两个称得到的结果之间进行估计,得到较为理想更科学合理的结果。这个过程我们也叫做数据融合。融合的数学表达式为:
z
^
=
K
z
1
+
(
1
−
K
)
z
2
\hat{z}=K z_{1}+(1-K) z_{2}
z^=Kz1+(1−K)z2
由于
z
1
z_{1}
z1和
z
2
z_{2}
z2相互独立,因此
z
^
\hat{z}
z^同样满足高斯分布,并且具有一下关系:
{
E
(
z
^
)
=
E
(
K
z
1
+
(
1
−
K
)
z
2
)
=
K
E
(
z
1
)
+
(
1
−
K
)
E
(
z
2
)
D
(
z
^
)
=
D
(
K
z
1
+
(
1
−
K
)
z
2
)
=
K
2
D
(
z
1
)
+
(
1
−
K
)
2
D
(
z
2
)
\left\{\begin{array}{l} E(\hat{z})=E\left(K z_{1}+(1-K) z_{2}\right)=K E\left(z_{1}\right)+(1-K) E\left(z_{2}\right) \\ D(\hat{z})=D\left(K z_{1}+(1-K) z_{2}\right)=K^{2} D\left(z_{1}\right)+(1-K)^{2} D\left(z_{2}\right) \end{array}\right.
{E(z^)=E(Kz1+(1−K)z2)=KE(z1)+(1−K)E(z2)D(z^)=D(Kz1+(1−K)z2)=K2D(z1)+(1−K)2D(z2)
显然,
K
K
K对于估计的结果有着关键的影响,因此我们的研究目的就转换成当
K
K
K取什么时,可以使得
z
^
\hat{z}
z^的方差
D
(
z
^
)
D(\hat{z})
D(z^)最小,并且求其对应的期望值
E
(
z
^
)
E(\hat{z})
E(z^),而最优估计
z
^
o
p
t
i
m
a
l
=
E
(
z
^
)
\hat{z}_{optimal}=E(\hat{z})
z^optimal=E(z^)。
为了实现上述研究目的,可以用导数的性质来求解
K
K
K的值,显然是关于
K
K
K的函数,其关系为:
D
(
z
^
)
=
K
2
D
(
z
1
)
+
(
1
−
K
)
2
D
(
z
2
)
D
(
z
^
)
=
0.2
K
2
+
(
1
−
2
K
+
K
2
)
×
0.4
D
(
z
^
)
=
0.6
K
2
−
0.8
K
+
0.4
D(\hat{z})=K^{2} D\left(z_{1}\right)+(1-K)^{2} D\left(z_{2}\right) \\ D(\hat{z})=0.2 K^{2}+\left(1-2 K+K^{2}\right) \times 0.4 \\ D(\hat{z})=0.6 K^{2}-0.8 K+0.4
D(z^)=K2D(z1)+(1−K)2D(z2)D(z^)=0.2K2+(1−2K+K2)×0.4D(z^)=0.6K2−0.8K+0.4
对K求导并求极值点
D
′
(
z
^
)
=
1.2
K
−
0.8
=
0
D^{\prime}(\hat{z})=1.2 K-0.8=0
D′(z^)=1.2K−0.8=0
可以解得:
K
=
2
3
K=\frac{2}{3}
K=32
综上,当
K
=
2
3
K=\frac{2}{3}
K=32时,得到的估计值分布的方差最小,此时期望值为
z
^
o
p
t
i
m
a
l
=
E
(
z
^
)
∣
K
=
2
3
=
5.1
\hat{z}_{optimal}=E(\hat{z})|_{K=\frac{2}{3}}=5.1
z^optimal=E(z^)∣K=32=5.1
当然了,原本也思考把和对调是否结果一致,即当估计算式转变为以下时,最优估计结果是否一致。
z
^
=
K
z
2
+
(
1
−
K
)
z
1
\hat{z}=K z_{2}+(1-K) z_{1}
z^=Kz2+(1−K)z1
通过上述的例子开始来讨论卡尔曼滤波器的内容。
卡尔曼滤波器
先给出卡尔曼滤波器的五条公式,即使不理解每条公式代表的具体数学含义,也希望能理解他的作用是什么,了解整个预测估计得思路是如何实现的。
Step1:用上一次的后验误差协方差矩阵计算先验误差协方差矩阵
P
k
+
1
−
=
A
P
k
A
T
+
Q
{\color{red} P_{k+1}^{-}=A P_{k} A^{T}+Q}
Pk+1−=APkAT+Q
Step2:先验估计
x
^
k
+
1
−
=
A
x
^
k
+
B
u
k
{\color{red} \hat{x}_{k+1}^{-}=A\hat{x}_{k}+Bu_{k}}
x^k+1−=Ax^k+Buk
Step3:计算卡尔曼增益
K
k
+
1
=
P
k
+
1
−
H
T
H
P
k
+
1
−
H
T
+
R
{\color{red} K_{k+1}=\frac{P_{k+1}^{-}H^{T}}{HP_{k+1}^{-}H^{T}+R}}
Kk+1=HPk+1−HT+RPk+1−HT
Step4:后验估计
x
^
k
+
1
=
x
^
k
+
1
−
+
K
k
+
1
(
z
k
+
1
−
H
x
^
k
+
1
−
)
{\color{red} \hat{x}_{k+1}=\hat{x}_{k+1}^{-}+K_{k+1}(z_{k+1}-H\hat{x}_{k+1}^{-})}
x^k+1=x^k+1−+Kk+1(zk+1−Hx^k+1−)
Step5:计算后验误差协方差矩阵
P
k
+
1
=
(
I
−
K
k
+
1
H
)
P
k
+
1
−
{\color{red} P_{k+1}=(I-K_{k+1}H)P_{k+1}^{-}}
Pk+1=(I−Kk+1H)Pk+1−
公式推导
对于一个离散系统中的状态向量,其状态空间转移方程为:
x
k
+
1
=
A
x
k
+
B
u
k
+
w
k
x_{k+1}=Ax_{k}+Bu_{k}+w_{k}
xk+1=Axk+Buk+wk
从传感器的角度来获得系统状态,其数学描述为:
z
k
+
1
=
H
x
k
+
1
+
v
k
+
1
z_{k+1}=Hx_{k+1}+v_{k+1}
zk+1=Hxk+1+vk+1
其中,
w
k
w_{k}
wk为
k
k
k时刻系统噪声,
v
k
v_{k}
vk为
k
k
k测量噪声,都满足高斯分布,即
w
k
∼
N
(
0
,
Q
)
v
k
∼
N
(
0
,
R
)
w_{k} \sim N(0,Q) \\ v_{k} \sim N(0,R)
wk∼N(0,Q)vk∼N(0,R)
其中,
Q
Q
Q和
R
R
R为两个分布的协方差矩阵。
无论是状态空间转移方程还是传感测测量转移方程,都包含了许多不确定性因素的干扰,这种干扰导致无法通过数学模型或者传感器获取到精准的状态向量,也就是说我们只能通过数学模型和传感器获取到两组不是很准确的状态向量。这种场景下就需要用到“数据融合”的知识。将这两组数据融合得到最为理想的估计值,也就是最优估计。
在对公式进行推导之前,我觉得有必要浅析一下对系统噪声和测量噪声的理解。假设系统十分理想,不存在噪声干扰。在这种情况下,我们可以直接了当的通过状态转移方程计算出下一时刻的状态,也可以通过传感器直接测量得到。这个时候无论是系统模型进行状态转移得到的状态向量还是通过传感器解算得到的状态向量,都是一个确定的,准确的量,也就不存在什么概率分布。
但是正因为现实生活中噪声无处不在,噪声是一种不稳定不确定的随机扰动。这种扰动是满足高斯分布的一种概率分布。因此对于不考虑噪声扰动的两个估计值——系统状态转移的状态向量估计值
x
^
k
+
1
−
\hat{x}_{k+1}^{-}
x^k+1−(也就是先验估计)和传感器解算的状态向量估计值
x
^
k
+
1
M
E
A
\hat{x}_{k+1}^{MEA}
x^k+1MEA,其表达式为:
x
^
k
+
1
−
=
A
x
^
k
+
B
u
k
x
^
k
+
1
M
E
A
=
H
−
z
k
+
1
{\color{blue} \hat{x}_{k+1}^{-}=A\hat{x}_{k}+Bu_{k}} \\ \hat{x}_{k+1}^{MEA}=H^{-}z_{k+1}
x^k+1−=Ax^k+Bukx^k+1MEA=H−zk+1
因此,对于两组估计值来说,都是不太准确的数据。当我们对这两组数据进行数据融合的时候,我们可以回忆前面讲述数据融合时的处理,作出最优估计表达式为:
x
^
k
+
1
=
x
^
k
+
1
−
+
G
(
x
^
k
+
1
M
E
A
−
x
^
k
+
1
−
)
\hat{x}_{k+1}=\hat{x}_{k+1}^{-}+G(\hat{x}_{k+1}^{MEA}-\hat{x}_{k+1}^{-})
x^k+1=x^k+1−+G(x^k+1MEA−x^k+1−)
代入式子
x
^
k
+
1
M
E
A
=
H
−
z
k
+
1
\hat{x}_{k+1}^{MEA}=H^{-}z_{k+1}
x^k+1MEA=H−zk+1,可以得到
x
^
k
+
1
=
x
^
k
+
1
−
+
G
(
H
−
z
k
+
1
−
x
^
k
+
1
−
)
\hat{x}_{k+1}=\hat{x}_{k+1}^{-}+G(H^{-}z_{k+1}-\hat{x}_{k+1}^{-})
x^k+1=x^k+1−+G(H−zk+1−x^k+1−)
由于公式中存在逆矩阵,为了简化运算,我们令
G
=
K
k
+
1
H
G=K_{k+1}H
G=Kk+1H,可以得到
x
^
k
+
1
=
x
^
k
+
1
−
+
K
k
+
1
H
(
H
−
z
k
+
1
−
x
^
k
+
1
−
)
=
>
x
^
k
+
1
=
x
^
k
+
1
−
+
K
k
+
1
(
z
k
+
1
−
H
x
^
k
+
1
−
)
\hat{x}_{k+1}=\hat{x}_{k+1}^{-}+K_{k+1}H(H^{-}z_{k+1}-\hat{x}_{k+1}^{-}) \\ =>{\color{blue} \hat{x}_{k+1}=\hat{x}_{k+1}^{-}+K_{k+1}(z_{k+1}-H\hat{x}_{k+1}^{-}) }
x^k+1=x^k+1−+Kk+1H(H−zk+1−x^k+1−)=>x^k+1=x^k+1−+Kk+1(zk+1−Hx^k+1−)
上面的公式就是后验估计的公式。下一步需要讨论一个问题:当
K
k
+
1
K_{k+1}
Kk+1取何值时,可以使得估计值达到最优。
我们需要思考几个问题:
①最优的标准是什么?
②这个标准的数学表达形式是什么?
在上面的讨论中,我们清楚了无论是先验估计
x
^
k
+
1
−
\hat{x}_{k+1}^{-}
x^k+1−还是传感器解算估计
x
^
k
+
1
M
E
A
\hat{x}_{k+1}^{MEA}
x^k+1MEA,都是满足一定条件的高斯分布,因此其加权和所得到的后验估计
x
^
k
+
1
\hat{x}_{k+1}
x^k+1,也是遵从一定条件的高斯分布。因此,我们可以通过调节
K
k
+
1
K_{k+1}
Kk+1的取值,来调整分布的特性,使得这个后验估计
x
^
k
+
1
\hat{x}_{k+1}
x^k+1的方差达到最小值,此时其数学期望就是最优的估计。
为了以数学形式更好的描述最优的标准,引入一个后验误差向量
e
k
+
1
e_{k+1}
ek+1,其表达式为:
e
k
+
1
=
x
k
+
1
−
x
^
k
+
1
e_{k+1}=x_{k+1}-\hat{x}_{k+1}
ek+1=xk+1−x^k+1
引入
x
^
k
+
1
=
x
^
k
+
1
−
+
K
k
+
1
(
z
k
+
1
−
H
x
^
k
+
1
−
)
\hat{x}_{k+1}=\hat{x}_{k+1}^{-}+K_{k+1}(z_{k+1}-H\hat{x}_{k+1}^{-})
x^k+1=x^k+1−+Kk+1(zk+1−Hx^k+1−)得到:
e
k
+
1
=
x
k
+
1
−
x
^
k
+
1
−
−
K
k
+
1
(
z
k
+
1
−
H
x
^
k
+
1
−
)
e_{k+1}=x_{k+1}-\hat{x}_{k+1}^{-}-K_{k+1}(z_{k+1}-H\hat{x}_{k+1}^{-})
ek+1=xk+1−x^k+1−−Kk+1(zk+1−Hx^k+1−)
引入
z
k
+
1
=
H
x
k
+
1
+
v
k
+
1
z_{k+1}=Hx_{k+1}+v_{k+1}
zk+1=Hxk+1+vk+1得到:
e
k
+
1
=
x
k
+
1
−
x
^
k
+
1
−
−
K
k
+
1
(
H
x
k
+
1
+
v
k
+
1
−
H
x
^
k
+
1
−
)
=
>
e
k
+
1
=
(
x
k
+
1
−
x
^
k
+
1
−
)
−
K
k
+
1
H
x
k
+
1
−
K
k
+
1
v
k
+
1
+
K
k
+
1
H
x
^
k
+
1
−
=
>
e
k
+
1
=
(
x
k
+
1
−
x
^
k
+
1
−
)
−
K
k
+
1
H
(
x
k
+
1
−
x
^
k
+
1
−
)
−
K
k
+
1
v
k
+
1
=
>
e
k
+
1
=
(
I
−
K
k
+
1
H
)
(
x
k
+
1
−
x
^
k
+
1
−
)
−
K
k
+
1
v
k
+
1
e_{k+1}=x_{k+1}-\hat{x}_{k+1}^{-}-K_{k+1}(Hx_{k+1}+v_{k+1}-H\hat{x}_{k+1}^{-}) \\ =>e_{k+1}=(x_{k+1}-\hat{x}_{k+1}^{-})-K_{k+1}Hx_{k+1}-K_{k+1}v_{k+1}+K_{k+1}H\hat{x}_{k+1}^{-} \\ =>e_{k+1}=(x_{k+1}-\hat{x}_{k+1}^{-})-K_{k+1}H(x_{k+1}-\hat{x}_{k+1}^{-})-K_{k+1}v_{k+1} \\ =>{\color{green} e_{k+1}=(I-K_{k+1}H)(x_{k+1}-\hat{x}_{k+1}^{-})-K_{k+1}v_{k+1}}
ek+1=xk+1−x^k+1−−Kk+1(Hxk+1+vk+1−Hx^k+1−)=>ek+1=(xk+1−x^k+1−)−Kk+1Hxk+1−Kk+1vk+1+Kk+1Hx^k+1−=>ek+1=(xk+1−x^k+1−)−Kk+1H(xk+1−x^k+1−)−Kk+1vk+1=>ek+1=(I−Kk+1H)(xk+1−x^k+1−)−Kk+1vk+1
显然,对于误差
e
k
+
1
e_{k+1}
ek+1来说,其表达式包含了先验估计
x
^
k
+
1
−
\hat{x}_{k+1}^{-}
x^k+1−以及测量噪声
v
k
+
1
v_{k+1}
vk+1两个高斯分布,因此
e
k
+
1
e_{k+1}
ek+1同样满足高斯分布。由于
x
^
k
+
1
−
\hat{x}_{k+1}^{-}
x^k+1−和
v
k
+
1
v_{k+1}
vk+1的数学期望均为0,因此
e
k
+
1
e_{k+1}
ek+1的数学期望也为0(这个地方可以用概率论数字期望中的定理简单验证得到),其协方差矩阵为
P
k
+
1
P_{k+1}
Pk+1。
因此,这个后验误差可以描述为:
e
k
+
1
∼
N
(
0
,
P
k
+
1
)
e_{k+1} \sim N(0,P_{k+1})
ek+1∼N(0,Pk+1)
那么,现在的问题是,只要尽可能的保证这个后验误差的协方差达到最小值,此时的卡尔曼增益
K
k
+
1
K_{k+1}
Kk+1就是最优解。那么对于后验误差协方差矩阵来说,其计算方式为:
P
k
+
1
=
E
(
e
k
+
1
e
k
+
1
T
)
P_{k+1}=E(e_{k+1}e_{k+1}^{T})
Pk+1=E(ek+1ek+1T)
插个眼,对于这个结论有兴趣了解推导过程的跳到下面协方差矩阵推导部分看。若已经清楚或不感兴趣可以忽略
令
e
k
+
1
−
=
x
k
+
1
−
x
^
k
+
1
−
e_{k+1}^{-}=x_{k+1}-\hat{x}_{k+1}^{-}
ek+1−=xk+1−x^k+1−,可以得到
e
k
+
1
=
(
I
−
K
k
+
1
H
)
e
k
+
1
−
−
K
k
+
1
v
k
+
1
e_{k+1}=(I-K_{k+1}H)e_{k+1}^{-}-K_{k+1}v_{k+1}
ek+1=(I−Kk+1H)ek+1−−Kk+1vk+1
那么可以得到
e
k
+
1
e
k
+
1
T
=
(
(
I
−
K
k
+
1
H
)
e
k
+
1
−
−
K
k
+
1
v
k
+
1
)
(
(
I
−
K
k
+
1
H
)
e
k
+
1
−
−
K
k
+
1
v
k
+
1
)
T
e_{k+1}e_{k+1}^{T}=((I-K_{k+1}H)e_{k+1}^{-}-K_{k+1}v_{k+1})((I-K_{k+1}H)e_{k+1}^{-}-K_{k+1}v_{k+1})^{T}
ek+1ek+1T=((I−Kk+1H)ek+1−−Kk+1vk+1)((I−Kk+1H)ek+1−−Kk+1vk+1)T
展开式子后可以得到
e
k
+
1
e
k
+
1
T
=
(
I
−
K
k
+
1
H
)
e
k
+
1
−
e
k
+
1
−
T
(
I
−
K
k
+
1
H
)
T
−
K
k
+
1
v
k
+
1
e
k
+
1
−
T
(
I
−
K
k
+
1
H
)
T
−
(
I
−
K
k
+
1
H
)
e
k
+
1
−
v
k
+
1
T
K
k
+
1
T
+
K
k
+
1
v
k
+
1
v
k
+
1
T
K
k
+
1
T
e_{k+1}e_{k+1}^{T}=(I-K_{k+1}H)e_{k+1}^{-}e_{k+1}^{-T}(I-K_{k+1}H)^{T} \\ -K_{k+1}v_{k+1}e_{k+1}^{-T}(I-K_{k+1}H)^{T}-(I-K_{k+1}H)e_{k+1}^{-}v_{k+1}^{T}K_{k+1}^{T} \\ +K_{k+1}v_{k+1}v_{k+1}^{T}K_{k+1}^{T}
ek+1ek+1T=(I−Kk+1H)ek+1−ek+1−T(I−Kk+1H)T−Kk+1vk+1ek+1−T(I−Kk+1H)T−(I−Kk+1H)ek+1−vk+1TKk+1T+Kk+1vk+1vk+1TKk+1T
令
P
k
+
1
−
=
E
(
e
k
+
1
−
e
k
+
1
−
T
)
P_{k+1}^{-}=E(e_{k+1}^{-}e_{k+1}^{-T})
Pk+1−=E(ek+1−ek+1−T),其中
P
k
+
1
−
P_{k+1}^{-}
Pk+1−为先验误差协方差,可以得到
E
(
e
k
+
1
e
k
+
1
T
)
=
(
I
−
K
k
+
1
H
)
E
(
e
k
+
1
−
e
k
+
1
−
T
)
(
I
−
K
k
+
1
H
)
T
−
K
k
+
1
E
(
v
k
+
1
e
k
+
1
−
T
)
(
I
−
K
k
+
1
H
)
T
−
(
I
−
K
k
+
1
H
)
E
(
e
k
+
1
−
v
k
+
1
T
)
K
k
+
1
T
+
K
k
+
1
E
(
v
k
+
1
v
k
+
1
T
)
K
k
+
1
T
E(e_{k+1}e_{k+1}^{T})=(I-K_{k+1}H)E(e_{k+1}^{-}e_{k+1}^{-T})(I-K_{k+1}H)^{T} \\ -K_{k+1}E(v_{k+1}e_{k+1}^{-T})(I-K_{k+1}H)^{T}-(I-K_{k+1}H)E(e_{k+1}^{-}v_{k+1}^{T})K_{k+1}^{T} \\ +K_{k+1}E(v_{k+1}v_{k+1}^{T})K_{k+1}^{T}
E(ek+1ek+1T)=(I−Kk+1H)E(ek+1−ek+1−T)(I−Kk+1H)T−Kk+1E(vk+1ek+1−T)(I−Kk+1H)T−(I−Kk+1H)E(ek+1−vk+1T)Kk+1T+Kk+1E(vk+1vk+1T)Kk+1T
其中,由于
v
k
+
1
v_{k+1}
vk+1和
e
k
+
1
e_{k+1}
ek+1相互独立,且各自的数学期望为0,则可以得到
E
(
v
k
+
1
e
k
+
1
)
=
E
(
v
k
+
1
)
E
(
e
k
+
1
)
=
0
E(v_{k+1}e_{k+1})=E(v_{k+1})E(e_{k+1})=0
E(vk+1ek+1)=E(vk+1)E(ek+1)=0
则原等式变换为:
E
(
e
k
+
1
e
k
+
1
T
)
=
(
I
−
K
k
+
1
H
)
P
k
+
1
−
(
I
−
K
k
+
1
H
)
T
+
K
k
+
1
R
K
k
+
1
T
E(e_{k+1}e_{k+1}^{T})=(I-K_{k+1}H)P_{k+1}^{-}(I-K_{k+1}H)^{T}+K_{k+1}RK_{k+1}^{T}
E(ek+1ek+1T)=(I−Kk+1H)Pk+1−(I−Kk+1H)T+Kk+1RKk+1T
再一次将式子展开可以得到
P
k
+
1
=
E
(
e
k
+
1
e
k
+
1
T
)
=
P
k
+
1
−
−
K
k
+
1
H
P
k
+
1
−
−
P
k
+
1
−
H
T
K
k
+
1
T
+
K
k
+
1
H
P
k
+
1
−
H
T
K
k
+
1
T
+
K
k
+
1
R
K
k
+
1
T
P_{k+1}=E(e_{k+1}e_{k+1}^{T})=P_{k+1}^{-}-K_{k+1}HP_{k+1}^{-}-P_{k+1}^{-}H^{T}K_{k+1}^{T}+K_{k+1}HP_{k+1}^{-}H^{T}K_{k+1}^{T}+K_{k+1}RK_{k+1}^{T}
Pk+1=E(ek+1ek+1T)=Pk+1−−Kk+1HPk+1−−Pk+1−HTKk+1T+Kk+1HPk+1−HTKk+1T+Kk+1RKk+1T
对于后验误差的协方差矩阵
P
k
+
1
P_{k+1}
Pk+1,我们并不关心其误差向量元素之间的协方差关系,我们只希望其每个元素的方差尽可能的小,对于一整个误差向量来说,其每个元素的方差之和达到最小时就是最优的状态。对于后验误差的协方差矩阵来说,向量各元素的方差之和就是该矩阵的迹,即:
∑
i
=
1
n
D
(
e
i
)
=
t
r
i
(
P
k
+
1
)
\sum_{i=1}^{n} D(e^{i})=tri(P_{k+1})
i=1∑nD(ei)=tri(Pk+1)
又该矩阵的迹与卡尔曼增益相关,因此也满足以下的关系:
t
r
i
(
P
k
+
1
)
=
f
(
K
k
+
1
)
tri(P_{k+1})=f(K_{k+1})
tri(Pk+1)=f(Kk+1)
对于找到最小值而言,我们只需要用迹的函数
t
r
i
(
P
k
+
1
)
tri(P_{k+1})
tri(Pk+1)对卡尔曼增益
K
k
+
1
K_{k+1}
Kk+1求导,即:
d
t
r
i
(
P
k
+
1
)
d
K
k
+
1
=
0
\frac{d tri(P_{k+1})}{dK_{k+1}}=0
dKk+1dtri(Pk+1)=0
注:这里为什么求导为0就可以定义为极小值而不是极大值,我是这么粗糙的理解的,从后验估计协方差矩阵的表达式:
P k + 1 = P k + 1 − − K k + 1 H P k + 1 − − P k + 1 − H T K k + 1 T + K k + 1 H P k + 1 − H T K k + 1 T + K k + 1 R K k + 1 T P_{k+1}=P_{k+1}^{-}-K_{k+1}HP_{k+1}^{-}-P_{k+1}^{-}H^{T}K_{k+1}^{T}+K_{k+1}HP_{k+1}^{-}H^{T}K_{k+1}^{T}+K_{k+1}RK_{k+1}^{T} Pk+1=Pk+1−−Kk+1HPk+1−−Pk+1−HTKk+1T+Kk+1HPk+1−HTKk+1T+Kk+1RKk+1T
可以看出,对于 K k + 1 K_{k+1} Kk+1来说,可以整理为:
P k + 1 = K k + 1 ( H P k + 1 − H T + R ) K k + 1 T − K k + 1 H P k + 1 − − P k + 1 − H T K k + 1 T + P k + 1 − P_{k+1}=K_{k+1}(HP_{k+1}^{-}H^{T}+R)K_{k+1}^{T}-K_{k+1}HP_{k+1}^{-}-P_{k+1}^{-}H^{T}K_{k+1}^{T}+P_{k+1}^{-} Pk+1=Kk+1(HPk+1−HT+R)Kk+1T−Kk+1HPk+1−−Pk+1−HTKk+1T+Pk+1−
可以等效的理解为:
P k + 1 = A K k + 1 2 − B K k + 1 + C P_{k+1}=AK_{k+1}^{2}-BK_{k+1}+C Pk+1=AKk+12−BKk+1+C
其中, A A A, B B B, C C C均为正数,因此其仅存在一个极小值。
注:这是我个人粗浅的理解,没有严格的数学证明,若有误请指正!!!
继续讨论这个公式
d
t
r
i
(
P
k
+
1
)
d
K
k
+
1
=
0
\frac{d tri(P_{k+1})}{dK_{k+1}}=0
dKk+1dtri(Pk+1)=0,对其展开后可以得到
d
t
r
i
(
P
k
+
1
)
d
K
k
+
1
=
d
t
r
i
(
P
k
+
1
−
)
d
K
k
+
1
−
d
t
r
i
(
K
k
+
1
H
P
k
+
1
−
)
d
K
k
+
1
−
d
t
r
i
(
P
k
+
1
−
H
T
K
k
+
1
T
)
d
K
k
+
1
+
d
t
r
i
(
K
k
+
1
H
P
k
+
1
−
H
T
K
k
+
1
T
)
d
K
k
+
1
+
d
t
r
i
(
K
k
+
1
R
K
k
+
1
T
)
d
K
k
+
1
⇒
d
t
r
i
(
P
k
+
1
)
d
K
k
+
1
=
d
t
r
i
(
P
k
+
1
−
)
d
K
k
+
1
−
2
d
t
r
i
(
K
k
+
1
H
P
k
+
1
−
)
d
K
k
+
1
+
d
t
r
i
(
K
k
+
1
H
P
k
+
1
−
H
T
K
k
+
1
T
)
d
K
k
+
1
+
d
t
r
i
(
K
k
+
1
R
K
k
+
1
T
)
d
K
k
+
1
\begin{array}{l} \frac{d t r i\left(P_{k+1}\right)}{d K_{k+1}}=\frac{d t r i\left(P_{k+1}^{-}\right)}{d K_{k+1}}-\frac{d t r i\left(K_{k+1} H P_{k+1}^{-}\right)}{d K_{k+1}}-\frac{d t r i\left(P_{k+1}^{-} H^{T} K_{k+1}^{T}\right)}{d K_{k+1}}+\frac{d t r i\left(K_{k+1} H P_{k+1}^{-} H^{T} K_{k+1}{ }^{T}\right)}{d K_{k+1}}+\frac{d t r i\left(K_{k+1} R K_{k+1}^{T}\right)}{d K_{k+1}} \Rightarrow \\ \frac{d t r i\left(P_{k+1}\right)}{d K_{k+1}}=\frac{d t r i\left(P_{k+1}^{-}\right)}{d K_{k+1}}-2 \frac{d t r i\left(K_{k+1} H P_{k+1}^{-}\right)}{d K_{k+1}}+\frac{d t r i\left(K_{k+1} H P_{k+1}^{-} H^{T} K_{k+1}{ }^{T}\right)}{d K_{k+1}}+\frac{d t r i\left(K_{k+1} R K_{k+1}^{T}\right)}{d K_{k+1}} \end{array}
dKk+1dtri(Pk+1)=dKk+1dtri(Pk+1−)−dKk+1dtri(Kk+1HPk+1−)−dKk+1dtri(Pk+1−HTKk+1T)+dKk+1dtri(Kk+1HPk+1−HTKk+1T)+dKk+1dtri(Kk+1RKk+1T)⇒dKk+1dtri(Pk+1)=dKk+1dtri(Pk+1−)−2dKk+1dtri(Kk+1HPk+1−)+dKk+1dtri(Kk+1HPk+1−HTKk+1T)+dKk+1dtri(Kk+1RKk+1T)
这里需要利用两个关于矩阵迹的结论:
d
t
r
i
(
A
B
)
d
A
=
B
T
d
tri
(
A
B
A
T
)
d
A
=
2
A
B
\begin{aligned} \frac{d t r i(A B)}{d A} &=B^{T} \\ \frac{d \operatorname{tri}\left(A B A^{T}\right)}{d A} &=2 A B \end{aligned}
dAdtri(AB)dAdtri(ABAT)=BT=2AB
代入后可以得到:
d
t
r
i
(
P
k
+
1
)
d
K
k
+
1
=
0
−
2
(
H
P
k
+
1
−
)
T
+
2
K
k
+
1
H
P
k
+
1
−
H
T
+
2
K
k
+
1
R
=
0
\frac{d t r i\left(P_{k+1}\right)}{d K_{k+1}}=0-2\left(H P_{k+1}^{-}\right)^{T}+2 K_{k+1} H P_{k+1}^{-} H^{T}+2 K_{k+1} R=0
dKk+1dtri(Pk+1)=0−2(HPk+1−)T+2Kk+1HPk+1−HT+2Kk+1R=0
解方程可以得
K
k
+
1
=
P
k
+
1
−
H
T
H
P
k
+
1
−
H
T
+
R
⇒
P
k
+
1
−
H
T
H
P
k
+
1
−
H
T
+
R
(
P
k
+
1
−
=
P
k
+
1
−
T
)
{\color{blue} K_{k+1}=\frac{P_{k+1}^{-} H^{T}}{H P_{k+1}^{-} H^{T}+R} \Rightarrow \frac{P_{k+1}^{-} H^{T}}{H P_{k+1}^{-} H^{T}+R}\left(P_{k+1}^{-}=P_{k+1}^{-}{ }^{T}\right)}
Kk+1=HPk+1−HT+RPk+1−HT⇒HPk+1−HT+RPk+1−HT(Pk+1−=Pk+1−T)
协方差矩阵是实对称矩阵
关于迹的定理公式可以直接通过矩阵运算的定义推导出来,此处不做推导。
从卡尔曼增益的定义来看,其中只有先验误差协方差矩阵
P
k
+
1
−
P_{k+1}^{-}
Pk+1−我们未求出,其他均为常矩阵可直接得到。接下来就开始计算
P
k
+
1
−
P_{k+1}^{-}
Pk+1−。
和后验误差协方差矩阵的定义同理,可以得到:
P
k
+
1
−
=
E
(
e
k
+
1
−
e
k
+
1
−
T
)
P_{k+1}^{-}=E(e_{k+1}^{-}e_{k+1}^{-T})
Pk+1−=E(ek+1−ek+1−T)
又由于
e
k
+
1
−
=
x
k
+
1
−
x
^
k
+
1
−
⇒
e
k
+
1
−
=
A
x
k
+
B
u
k
+
w
k
−
A
x
^
k
−
B
u
k
⇒
e
k
+
1
−
=
A
(
x
k
−
x
^
k
)
+
w
k
⇒
e
k
+
1
−
=
A
e
k
+
w
k
\begin{array}{l} e_{k+1}^{-}=x_{k+1}-\hat{x}_{k+1}^{-} \Rightarrow \\ e_{k+1}^{-}=A x_{k}+B u_{k}+w_{k}-A \hat{x}_{k}-B u_{k} \Rightarrow \\ e_{k+1}^{-}=A\left(x_{k}-\hat{x}_{k}\right)+w_{k} \Rightarrow \\ e_{k+1}^{-}=A e_{k}+w_{k} \end{array}
ek+1−=xk+1−x^k+1−⇒ek+1−=Axk+Buk+wk−Ax^k−Buk⇒ek+1−=A(xk−x^k)+wk⇒ek+1−=Aek+wk
代入后可以得到:
P
k
+
1
−
=
E
(
e
k
+
1
−
e
k
+
1
−
)
⇒
P
k
+
1
−
=
E
(
(
A
e
k
+
w
k
)
(
A
e
k
+
w
k
)
T
)
⇒
P
k
+
1
−
=
E
(
(
A
e
k
+
w
k
)
(
e
k
T
A
T
+
w
k
T
)
)
⇒
P
k
+
1
−
=
E
(
A
e
k
e
k
T
A
T
+
w
k
e
k
T
A
T
+
A
e
k
w
k
T
+
w
k
w
k
T
)
⇒
P
k
+
1
−
=
A
E
(
e
k
e
k
T
)
A
T
+
E
(
w
k
e
k
T
)
A
T
+
A
E
(
e
k
w
k
T
)
+
E
(
w
k
w
k
T
)
⇒
P
k
+
1
−
=
A
P
k
A
T
+
E
(
w
k
)
E
(
e
k
T
)
A
T
+
A
E
(
e
k
)
E
(
w
k
T
)
+
Q
⇒
P
k
+
1
−
=
A
P
k
A
T
+
Q
\begin{array}{l} P_{k+1}^{-}=E\left(e_{k+1}^{-} e_{k+1}^{-}\right) \Rightarrow \\ P_{k+1}^{-}=E\left(\left(A e_{k}+w_{k}\right)\left(A e_{k}+w_{k}\right)^{T}\right) \Rightarrow \\ P_{k+1}^{-}=E\left(\left(A e_{k}+w_{k}\right)\left(e_{k}^{T} A^{T}+w_{k}^{T}\right)\right) \Rightarrow \\ P_{k+1}^{-}=E\left(A e_{k} e_{k}^{T} A^{T}+w_{k} e_{k}^{T} A^{T}+A e_{k} w_{k}^{T}+w_{k} w_{k}^{T}\right) \Rightarrow \\ P_{k+1}^{-}=A E\left(e_{k} e_{k}^{T}\right) A^{T}+E\left(w_{k} e_{k}^{T}\right) A^{T}+A E\left(e_{k} w_{k}^{T}\right)+E\left(w_{k} w_{k}^{T}\right) \Rightarrow \\ P_{k+1}^{-}=A P_{k} A^{T}+E\left(w_{k}\right) E\left(e_{k}^{T}\right) A^{T}+A E\left(e_{k}\right) E\left(w_{k}^{T}\right)+Q \Rightarrow \\ {\color{blue} P_{k+1}^{-}=A P_{k} A^{T}+Q} \end{array}
Pk+1−=E(ek+1−ek+1−)⇒Pk+1−=E((Aek+wk)(Aek+wk)T)⇒Pk+1−=E((Aek+wk)(ekTAT+wkT))⇒Pk+1−=E(AekekTAT+wkekTAT+AekwkT+wkwkT)⇒Pk+1−=AE(ekekT)AT+E(wkekT)AT+AE(ekwkT)+E(wkwkT)⇒Pk+1−=APkAT+E(wk)E(ekT)AT+AE(ek)E(wkT)+Q⇒Pk+1−=APkAT+Q
从公式里可以看出,需要用到上一时刻的后验误差协方差矩阵,因此我们也需要进行后验误差协方差矩阵
P
k
+
1
P_{k+1}
Pk+1的推导。
在计算卡尔曼增益时已经推导过后验误差协方差矩阵的公式,即:
P
k
+
1
=
P
k
+
1
−
−
K
k
+
1
H
P
k
+
1
−
−
P
k
+
1
−
H
T
K
k
+
1
T
+
K
k
+
1
H
P
k
+
1
−
H
T
K
k
+
1
T
+
K
k
+
1
R
K
k
+
1
T
P_{k+1}=P_{k+1}^{-}-K_{k+1}HP_{k+1}^{-}-P_{k+1}^{-}H^{T}K_{k+1}^{T}+K_{k+1}HP_{k+1}^{-}H^{T}K_{k+1}^{T}+K_{k+1}RK_{k+1}^{T}
Pk+1=Pk+1−−Kk+1HPk+1−−Pk+1−HTKk+1T+Kk+1HPk+1−HTKk+1T+Kk+1RKk+1T
代入先验误差协方差矩阵
P
k
+
1
−
=
A
P
k
A
T
+
Q
P_{k+1}^{-}=A P_{k} A^{T}+Q
Pk+1−=APkAT+Q以及卡尔曼增益
K
k
+
1
=
P
k
+
1
−
H
T
H
P
k
+
1
−
H
T
+
R
K_{k+1}=\frac{P_{k+1}^{-} H^{T}}{H P_{k+1}^{-} H^{T}+R}
Kk+1=HPk+1−HT+RPk+1−HT,可以得到
P
k
+
1
=
P
k
+
1
−
−
K
k
+
1
H
P
k
+
1
−
−
P
k
+
1
−
H
T
K
k
+
1
T
+
P
k
+
1
−
H
T
H
P
k
+
1
−
H
T
+
R
(
H
P
k
+
1
−
H
T
+
R
)
K
k
+
1
T
⇒
P
k
+
1
=
P
k
+
1
−
−
K
k
+
1
H
P
k
+
1
−
−
P
k
+
1
−
H
T
K
k
+
1
T
+
P
k
+
1
−
H
T
K
k
+
1
T
⇒
P
k
+
1
=
P
k
+
1
−
−
K
k
+
1
H
P
k
+
1
−
⇒
P
k
+
1
=
(
I
−
K
k
+
1
H
)
P
k
+
1
−
\begin{array}{l} P_{k+1}=P_{k+1}^{-}-K_{k+1} H P_{k+1}^{-}-P_{k+1}^{-} H^{T} K_{k+1}{ }^{T}+\frac{P_{k+1}^{-} H^{T}}{H P_{k+1}^{-} H^{T}+R}\left(H P_{k+1}^{-} H^{T}+R\right) K_{k+1}{ }^{T} \Rightarrow \\ P_{k+1}=P_{k+1}^{-}-K_{k+1} H P_{k+1}^{-}-P_{k+1}^{-} H^{T} K_{k+1}{ }^{T}+P_{k+1}^{-} H^{T} K_{k+1}{ }^{T} \Rightarrow \\ P_{k+1}=P_{k+1}^{-}-K_{k+1} H P_{k+1}^{-} \Rightarrow \\ {\color{blue} P_{k+1}=\left(I-K_{k+1} H\right) P_{k+1}^{-}} \end{array}
Pk+1=Pk+1−−Kk+1HPk+1−−Pk+1−HTKk+1T+HPk+1−HT+RPk+1−HT(HPk+1−HT+R)Kk+1T⇒Pk+1=Pk+1−−Kk+1HPk+1−−Pk+1−HTKk+1T+Pk+1−HTKk+1T⇒Pk+1=Pk+1−−Kk+1HPk+1−⇒Pk+1=(I−Kk+1H)Pk+1−
从以上的推导中,我们整理一下卡尔曼滤波器的思路,主要分为预测和校正的两个部分。首先利用上一次的后验估计值和输入向量预测出当前时刻的先验估计,并计算出先验误差协方差矩阵。然后根据先验误差协方差矩阵计算卡尔曼增益,后验估计,并更新后验误差协方差矩阵。
预测
x ^ k + 1 − = A x ^ k + B u k \hat{x}_{k+1}^{-}=A\hat{x}_{k}+Bu_{k} x^k+1−=Ax^k+Buk
P k + 1 − = A P k A T + Q P_{k+1}^{-}=A P_{k} A^{T}+Q Pk+1−=APkAT+Q
第一步根据系统状态转移(个人理解为系统自身存在的惯性)和输入的影响来估计出系统下一时刻的状态
第二步计算出这一估计的误差协方差矩阵(个人理解为对于这一估计的信任程度,受到上一时刻后验估计可信程度和系统噪声影响)
校正
K k + 1 = P k + 1 − H T H P k + 1 − H T + R K_{k+1}=\frac{P_{k+1}^{-} H^{T}}{H P_{k+1}^{-} H^{T}+R} Kk+1=HPk+1−HT+RPk+1−HT
x ^ k + 1 = x ^ k + 1 − + K k + 1 ( z k + 1 − H x ^ k + 1 − ) \hat{x}_{k+1}=\hat{x}_{k+1}^{-}+K_{k+1}(z_{k+1}-H\hat{x}_{k+1}^{-}) x^k+1=x^k+1−+Kk+1(zk+1−Hx^k+1−)
P k + 1 = ( I − K k + 1 H ) P k + 1 − P_{k+1}=\left(I-K_{k+1} H\right) P_{k+1}^{-} Pk+1=(I−Kk+1H)Pk+1−
第一步是根据先验误差协方差矩阵计算卡尔曼增益系数(个人理解这一过程是在权衡先验估计和测算估计之间的比重,受到先验误差协方差矩阵和测量噪声的影响)
第二步是计算后验估计
第三步是计算出本次后验估计的误差协方差矩阵,为下一次预测做准备。
从公式上看,卡尔曼滤波属于递归性质的算法,因此对具体的模型预测时,需要设定好初始状态向量以及对应的后验误差协方差矩阵的初始值。这里初始状态向量容易理解,即初始时刻的系统状态,那此时的后验估计协方差矩阵意味着对这个初始状态的可靠性程度。因为初始状态向量也不是绝对准确的,存在一定偏差的,那么这个初始的后验误差协方差矩阵就代表这个意思。但在实际的系统中,后验误差协方差矩阵初始值无论多少,最终都会依据系统噪声和测量噪声进行收敛。
卡尔曼滤波器是建立在线性离散系统的基础上的,若系统为非线性系统,则不能直接套用以上结论,需要采用扩展卡尔曼滤波器,将系统线性化后再进行预测校正,其整体的思路是一致的。
协方差矩阵推导
对于一个变量
X
X
X来说,其方差表示为:
D
(
X
)
=
E
(
X
2
)
−
E
2
(
X
)
D(X)=E(X^{2})-E^{2}(X)
D(X)=E(X2)−E2(X)
对于两个变量
X
X
X
Y
Y
Y来说,其协方差表示为:
C
o
v
(
X
,
Y
)
=
E
(
X
Y
)
−
E
(
X
)
E
(
Y
)
Cov(X,Y)=E(XY)-E(X)E(Y)
Cov(X,Y)=E(XY)−E(X)E(Y)
从表达式来看,无论是方差还是协方差,最后都可以转化成数学期望的计算。那么对于一个向量
z
=
[
z
1
z
2
⋮
z
n
]
z=\begin{bmatrix} z_{1}\\z_{2} \\\vdots \\z_{n} \end{bmatrix}
z=⎣
⎡z1z2⋮zn⎦
⎤,其协方差矩阵
P
P
P应该表示为:
P
=
[
D
(
z
1
)
C
o
v
(
z
1
,
z
2
)
⋯
C
o
v
(
z
1
,
z
n
)
C
o
v
(
z
2
,
z
1
)
D
(
z
2
)
⋯
C
o
v
(
z
2
,
z
n
)
⋮
⋮
⋱
⋮
C
o
v
(
z
n
,
z
1
)
⋮
⋯
D
(
z
n
)
]
P=\begin{bmatrix} D(z_{1})& Cov(z_{1},z_{2})& \cdots & Cov(z_{1},z_{n})\\ Cov(z_{2},z_{1})& D(z_{2})& \cdots& Cov(z_{2},z_{n})\\ \vdots & \vdots& \ddots& \vdots\\ Cov(z_{n},z_{1})& \vdots& \cdots & D(z_{n}) \end{bmatrix}
P=⎣
⎡D(z1)Cov(z2,z1)⋮Cov(zn,z1)Cov(z1,z2)D(z2)⋮⋮⋯⋯⋱⋯Cov(z1,zn)Cov(z2,zn)⋮D(zn)⎦
⎤
代入
D
(
X
)
=
E
(
X
2
)
−
E
2
(
X
)
D(X)=E(X^{2})-E^{2}(X)
D(X)=E(X2)−E2(X)和
C
o
v
(
X
,
Y
)
=
E
(
X
Y
)
−
E
(
X
)
E
(
Y
)
Cov(X,Y)=E(XY)-E(X)E(Y)
Cov(X,Y)=E(XY)−E(X)E(Y),可以得到
P
=
[
E
(
z
1
2
)
−
E
2
(
z
1
)
E
(
z
1
z
2
)
−
E
(
z
1
)
E
(
z
2
)
⋯
E
(
z
1
z
n
)
−
E
(
z
1
)
E
(
z
n
)
E
(
z
2
z
1
)
−
E
(
z
2
)
E
(
z
1
)
E
(
z
2
2
)
−
E
2
(
z
2
)
⋯
E
(
z
2
z
n
)
−
E
(
z
2
)
E
(
z
n
)
⋮
⋮
⋱
⋮
E
(
z
n
z
1
)
−
E
(
z
n
)
E
(
z
1
)
⋮
⋯
E
(
z
n
2
)
−
E
2
(
z
n
)
]
P=\begin{bmatrix} E(z_{1}^{2})-E^{2}(z_{1})& E(z_{1}z_{2})-E(z_{1})E(z_{2})& \cdots & E(z_{1}z_{n})-E(z_{1})E(z_{n})\\ E(z_{2}z_{1})-E(z_{2})E(z_{1})& E(z_{2}^{2})-E^{2}(z_{2})& \cdots& E(z_{2}z_{n})-E(z_{2})E(z_{n})\\ \vdots & \vdots& \ddots& \vdots\\ E(z_{n}z_{1})-E(z_{n})E(z_{1})& \vdots& \cdots & E(z_{n}^{2})-E^{2}(z_{n}) \end{bmatrix}
P=⎣
⎡E(z12)−E2(z1)E(z2z1)−E(z2)E(z1)⋮E(znz1)−E(zn)E(z1)E(z1z2)−E(z1)E(z2)E(z22)−E2(z2)⋮⋮⋯⋯⋱⋯E(z1zn)−E(z1)E(zn)E(z2zn)−E(z2)E(zn)⋮E(zn2)−E2(zn)⎦
⎤
当向量
z
z
z的数学期望均为0时,则原式可以转化为:(先验误差和后验误差的数学期望均为0)
P
=
[
E
(
z
1
2
)
E
(
z
1
z
2
)
⋯
E
(
z
1
z
n
)
E
(
z
2
z
1
)
E
(
z
2
2
)
⋯
E
(
z
2
z
n
)
⋮
⋮
⋱
⋮
E
(
z
n
z
1
)
⋮
⋯
E
(
z
n
2
)
]
=
E
[
z
1
2
z
1
z
2
⋯
z
1
z
n
z
2
z
1
z
2
2
⋯
z
2
z
n
⋮
⋮
⋱
⋮
z
n
z
1
⋮
⋯
z
n
2
]
=
E
(
z
z
T
)
P=\begin{bmatrix} E(z_{1}^{2})& E(z_{1}z_{2})& \cdots & E(z_{1}z_{n})\\ E(z_{2}z_{1})& E(z_{2}^{2})& \cdots& E(z_{2}z_{n})\\ \vdots & \vdots& \ddots& \vdots\\ E(z_{n}z_{1})& \vdots& \cdots & E(z_{n}^{2}) \end{bmatrix} =E\begin{bmatrix} z_{1}^{2}& z_{1}z_{2}& \cdots & z_{1}z_{n}\\ z_{2}z_{1}& z_{2}^{2}& \cdots& z_{2}z_{n}\\ \vdots & \vdots& \ddots& \vdots\\ z_{n}z_{1}& \vdots& \cdots & z_{n}^{2} \end{bmatrix}=E(zz^{T})
P=⎣
⎡E(z12)E(z2z1)⋮E(znz1)E(z1z2)E(z22)⋮⋮⋯⋯⋱⋯E(z1zn)E(z2zn)⋮E(zn2)⎦
⎤=E⎣
⎡z12z2z1⋮znz1z1z2z22⋮⋮⋯⋯⋱⋯z1znz2zn⋮zn2⎦
⎤=E(zzT)
总结
借用一句话来说,数学是没办法科普的。所以尽管参杂了很多浅显易懂的文字表达,依旧不能完全表达出我对于卡尔曼滤波器的理解,更多是一种思考问题解决问题的逻辑。后续计划写一篇卡尔曼滤波器在实际工程中的设计和应用过程,以此来更加直观的体现出卡尔曼滤波器所发挥的作用。
在此鸣谢B站up主 DR_CAN 的相关视频