什么是卡尔曼滤波器#
卡尔曼滤波是一种递归的估计,即只要获知上一时刻状态的估计值以及当前状态的观测值就可以计算出当前状态的估计值,因此不需要记录观测或者估计的历史信息。
卡尔曼增益 K 是把预测/先验和观测/测量按不确定性加权融合的矩阵权重:它根据预测的不确定性和测量噪声决定“该信任多少观测、该信任多少预测”,从而把测量残差(innovation)映射回状态空间来修正预测。
卡尔曼滤波器的数学模型#
设离散线性系统模型有:
状态传播模型:观测模型:xk=Fkxk−1+Bkuk+wk,wk∼N(0,Qk)zk=Hkxk+vk,vk∼N(0,Rk)
- 状态传播模型是客观存在的,永远准确的,但无法直接观测的。
- 观测模型是通过测量手段得到的测量数据,包含噪声。
先验估计(状态预测):x^k∣k−1=Fkx^k−1∣k−1+Bkuk
先验协方差矩阵(先验估计的精度描述):Pk∣k−1=FkPk−1∣k−1FkT+Qk
令:Sk=HkPk∣k−1Hk⊤+Rk
卡尔曼增益:
Kk=Pk∣k−1Hk⊤Sk−1
后验估计(对先验估计进行更新):
x^k∣k=x^k∣k−1+Kkyk,yk=zk−Hkx^k∣k−1
后验协方差矩阵(后验估计的精度描述)Pk∣k=(I−KkHk)Pk∣k−1
很好,这里五个公式里面有10个神秘字母,下面将一一介绍他们各自的含义。
符号说明#
假设系统状态空间的维度为 n,测量(观测)空间的维度为 m
TIP所有带 ^ 的值都是估计值,不带 ^ 的都是真实值
状态空间#
真实状态向量#
xk:系统在时刻 k 的真实状态,无法观测(受噪声影响)维度 n×1。
预测状态/先验估计#
x^k∣k−1:在 k 时刻使用 k−1 时刻的后验估计和系统的状态转移矩阵得到的预测,叫先验估计,维度 n×1。
更新状态/后验估计#
x^k∣k:在包含 k 时刻观测后测量值的估计后验状态估计,维度 n×1。
测量(观测)空间#
测量噪声#
vk:m×1。测量噪声假设服从均值为0,协方差为 Rk 的高斯分布。
观测矩阵#
Hk:把状态(x^和x)映射到测量空间(z),m×n。观测矩阵每一行代表一个测量值,每一列代表一个状态值,矩阵中的元素代表状态对测量的影响程度。
理想测量向量#
zk:m×1。观测矩阵作用于真实状态向量
估计测量向量#
z^k:m×1。观测矩阵作用于先验估计(预测状态)
估计状态向量#
x^k:n×1。由测量值映射回状态空间后的估计状态向量。
测量残差#
yk:理想测量和估计测量间的误差。
yk=zk−z^k=zk−Hkx^k∣k−1,大小 m×1。 其中 zk 是理想测量向量。
-
(理想的和估计的)测量向量的维度取决于观测矩阵的行数,即 m。
-
测量是存在误差的,有噪声
-
观测矩阵作用于(xk)并加上测量噪声(vk),得到理想的测量空间(zk)
-
观测矩阵作用于(x^k)并加上测量噪声(vk),得到估计的测量空间(z^k)
-
理想测量向量 zk 和估计测量向量 z^k 由测量残差 yk 关联
NOTE比如说对于一个运动的物体,有速度,位置,加速度三个状态,则 k 时刻的 真实状态空间 xk(3×1,n=3)=pva ,但是观测手段只允许测量位置和速度,那么观测矩阵 Hk(2×3,m=2,n=3)=[100100] ,则 k 时刻得到的 理想测量空间 zk(m=2,n=1):
zk=Hk⋅xk+vk=[100100]⋅pva+vk=[pv]+vk
两个协方差矩阵#
WARNING不要尝试用严格的协方差矩阵数学定义去看待卡尔曼滤波中所谓的协方差矩阵
概率论中的协方差矩阵:
Σ=Cov(X)=E[(X−μ)(X−μ)T]X减的是均值,描述的是随机向量相对于它的真实均值的离散程度和分量间相关性。
卡尔曼滤波中的协方差矩阵:(这里以先验协方差为例)
Pk∣k−1=E[(xk−x^k∣k−1)(xk−x^k∣k−1)T]减的是估计值,描述的是真实状态xk(随机向量)相对于它的先验估计(估计值)xk∣k−1的估计精度(离散程度)和分量间相关性。
TIP简而言之,卡尔曼滤波中的协方差矩阵就是估计精度。
- 为啥不用真实的均值呢?因为这玩意没法观测。
- 这种做法是合理的,卡尔曼滤波器在线性、高斯、无偏初始条件下,保证了E[xk]=x^k∣k(证明过程略,可以参考维基),所以用估计值代替真实均值是合理的。
ek∣k−1=xk−x^k∣k−1:先验估计误差:真实状态与先验估计的误差。
ek∣k=xk−x^k∣k:后验误差估计:真实状态与后验估计的误差。
Pk∣k−1:先验协方差矩阵(预测精度):=E[ek∣k−1ek∣k−1T],n×n。
Pk∣k:后验协方差矩阵(更新精度):=E[ek∣kek∣kT],n×n。
从上文可知,估计分为先验估计 x^k∣k−1 和后验估计 x^k∣k ,这两个估计的”真实值”都是真实状态向量 xk 。Pk∣k−1 先验协方差矩阵描述先验估计的精度(预测的精度),Pk∣k 后验协方差矩阵描述后验估计的精度(更新的精度)。
两个协方差矩阵的推导
Pk∣k−1Pk∣k=E[ek∣k−1ek∣k−1T]=E[ek∣kek∣kT]1. 先验误差协方差 Pk∣k−1#
xkx^k∣k−1ek∣k−1Pk∣k−1∴Pk∣k−1系统的状态转移方程(描述客观的状态转移):=Fkxk−1+Bkuk+wk先验估计方程:=Fkx^k−1∣k−1+Bkuk(2)−(1)=xk−x^k∣k−1=(Fkxk−1+Bkuk+wk)−(Fkx^k−1∣k−1+Bkuk)=Fkek−1∣k−1+wk=E[ek∣k−1ek∣k−1⊤]=E[Fkek−1∣k−1ek−1∣k−1⊤Fk⊤]+E[Fkek−1∣k−1wk⊤]+E[wkek−1∣k−1⊤Fk⊤]E[wkwk⊤]因为 wk 独立,且零均值:E(wk)=0=FkE[ek−1∣k−1ek−1∣k−1⊤]Fk⊤+E[wkwk⊤]=FkPk−1∣k−1Fk⊤+Qk这里 Qk 是过程噪声的协方差。⋯(1)⋯(2)2. 后验误差协方差 Pk∣k#
x^k∣kzkek∣kPk∣k后验估计方程:=x^k∣k−1+Kk(zk−Hkx^k∣k−1)⋯(3)观测方程:=Hkxk+vk⋯(4)后验误差:=xk−x^k∣k=(4)→(3)=xk−x^k∣k−1−Kk(Hkxk+vk−Hkx^k∣k−1)=(I−KkHk)(xk−x^k∣k−1)−Kkvk=(I−KkHk)ek∣k−1−Kkvk=E[ek∣kek∣k⊤]=E[((I−KkHk)ek∣k−1−Kkvk)((I−KkHk)ek∣k−1−Kkvk)⊤]=(I−KkHk)E[ek∣k−1ek∣k−1⊤](I−KkHk)⊤+KkE[vkvk⊤]Kk⊤上一步用了协方差矩阵的简化算法的思路,即Σ=E[XXT]−μμT=(I−KkHk)Pk∣k−1(I−KkHk)⊤+KkRkKk⊤这里 Rk 是测量噪声的协方差矩阵(可选)化简到常见简化式的完整代数步骤#
常见简化形式为
Pk∣k=(I−KkHk)Pk∣k−1.下面给出它从 Joseph 形式代数化简的步骤,并说明何时成立。
从 Joseph 形式展开:
Pk∣k=Pk∣k−1−KkHkPk∣k−1−Pk∣k−1Hk⊤Kk⊤+KkHkPk∣k−1Hk⊤Kk⊤+KkRkKk⊤.把最后两项合并:
KkHkPk∣k−1Hk⊤Kk⊤+KkRkKk⊤=Kk(HkPk∣k−1Hk⊤+Rk)Kk⊤.记
Sk=HkPk∣k−1Hk⊤+Rk.若使用标准卡尔曼增益
Kk=Pk∣k−1Hk⊤Sk−1,则
KkSkKk⊤=(Pk∣k−1Hk⊤Sk−1)Sk(Sk−1HkPk∣k−1)=Pk∣k−1Hk⊤Sk−1HkPk∣k−1.注意到
Pk∣k−1Hk⊤Sk−1=Kk,且Sk−1HkPk∣k−1=Kk⊤.于是合并项退化为:
KkSkKk⊤=Pk∣k−1Hk⊤Kk⊤.回代 Joseph 展开式:
Pk∣k=Pk∣k−1−KkHkPk∣k−1−Pk∣k−1Hk⊤Kk⊤+Pk∣k−1Hk⊤Kk⊤=Pk∣k−1−KkHkPk∣k−1=(I−KkHk)Pk∣k−1.因此在采用标准卡尔曼增益并利用定义 Sk 后,Joseph 形式代数上会化简为常用的简化形式。但需注意:
- 简化式 Pk∣k=(I−KkHk)Pk∣k−1 并不显式包含 Rk,其成立依赖于上面用到的 Kk 的表达式(即利用了 KkSkKk⊤ 的代数关系)。
- 为了数值对称性与稳定性,工程上更常用 Joseph 形式带两个 A 的对称形式 KRK⊤,因为它显式保证了 Pk∣k 的对称正半定性。
这说明后验协方差是先验协方差经过“卡尔曼增益修正”后的结果。
上面说卡尔曼矩阵描述了估计值的估计精度,所以也能决定卡尔曼增益 Kk 的大小。这个在卡尔曼增益的解析中给出推导
真实状态空间中的状态转移相关的矩阵#
状态转移矩阵#
Fk:系统动力学,n×n。
这玩意描述真实系统的状态转移方式,即 xk−1⟶Fkxk。受其直接影响的量是真实状态 xk 和先验协方差 Pk∣k−1。
控制输入矩阵&控制向量#
Bk:n×p。
uk:p×1。
Bk 和 uk 描述系统受到的外部控制输入对状态的影响,比如汽车的油门、刹车、方向盘等。
过程噪声协方差#
Qk:n×n。
过程噪声协方差矩阵 Qk 表示系统在状态转移过程中可能存在的随机扰动或不确定性。比如风吹、路面不平等因素都会影响物体的运动状态,这个在 k 时刻对物体的不确定影响我们就用 Qk 来表示
测量噪声协方差#
Rk:m×m。
测量噪声协方差矩阵 Rk 表示测量仪器本身存在的误差或不确定性,比如测量仪器的精度有限,可能会有一定的误差,这个误差我们就用 Rk 来表示。
创新(残差)协方差#
Sk:也叫测量预测协方差。m×m,它表示把预测的不确定性投影到测量空间后,与测量噪声合并得到的“测量不确定性”。
卡尔曼增益的直观意义#
Kk:卡尔曼增益,大小为 n×m
用现代汉语描述一下卡尔曼增益在后验估计方程(状态更新)中的作用:
当前的估计值(后验估计)x^k∣k=上一次的估计(先验估计)x^k∣k−1+系数(卡尔曼增益)Kk⋅当前测量值-上一次估计值(zk−Hk⋅x^k∣k−1)卡尔曼滤波认为,后验估计 x^k∣k 是先验估计和测量值的线性组合,也就是进行数据融合,卡尔曼增益决定了”相信多少先验估计,相信多少测量值”
如果纯用数据融合的方式,那么公式的形式应该长这样
x^k∣k=x^k∣k−1+G⋅(x^k−x^k∣k−1)其中x^k=Hk−zkx^k 是根据测量值映射回状态空间后对客观状态的估计
G 是数据融合的系数矩阵,范围是 [0,I],当 G=0 时,完全相信先验估计 x^k∣k−1;当 G=I 时,完全相信测量值 x^k。
对系数 G 稍微变形一下就得到了标准的卡尔曼滤波公式:
令G=KkHk⇒x^k∣k=x^k∣k−1+Kk⋅(zk−Hk⋅x^k∣k−1)这里的 Kk 就是卡尔曼增益。范围是 [0,Hk−1],当 Kk=0 时,完全相信先验估计 x^k∣k−1;当 Kk=Hk−1 时,完全相信测量值 x^k。
下面推导可能有误,谨慎参考
卡尔曼增益的推导
这一坨东西的推导和后验协方差矩阵息息相关,后验协方差矩阵越“小”越好。对于这个“小”的标准,最常见的是后验协方差矩阵的迹(trace),也就是对角线各元素之和(对角线元素是后验误差的方差)
具体推导里面有矩阵求导bulabula之类的,主包的线性代数功底不好。
为了简化,记
P≡Pk∣k−1,K≡Kk,H=Hk,R=Rk结合之前后验协方差矩阵的推导过程,Joshep展开:
Pk∣k=(I−KH)P(I−KH)⊤+KRK⊤=(I−KH)P(I⊤−H⊤K⊤)+KRK⊤=P−KHP−PH⊤K⊤+KHPH⊤K⊤+KRK⊤其中 P 是已知的先验协方差矩阵,H 是已知的观测矩阵,R 是已知的测量噪声协方差矩阵,未知的是 K。
以迹为代价函数,对角线元素之和越小,说明后验误差的方差越小,估计精度越高
J(K)=tr(Pk∣k).用迹的性质展开,关于迹的性质,可以参考这篇博客。
迹的循环不变性 tr(ABC)=tr(BCA)
转置不改变迹 tr(P⊤)=tr(P),当然,P 本来就是对称的。
J(K)=tr((I−KH)P(I−KH)⊤)+tr(KRK⊤)=tr(P)−2tr(KHP)+tr(KHPH⊤K⊤)+tr(KRK⊤).为了处理后面两项,定义
S=HPH⊤+R,则
J(K)=tr(P)−2tr(KHP)+tr(KSK⊤)矩阵微分与梯度
对 K 求导。用到两个常用公式:
-
若 A 为常矩阵∂K∂tr(KA)=A⊤。
-
若 A 对称,∂K∂tr(KAK⊤)=2KA。
于是
∂K∂J(K)=−2(HP)⊤+2KS=−2PH⊤+2KS.optimality 条件与最优解,令梯度为零:
−2PH⊤+2KS=0⟹KS=PH⊤.K=PH⊤S−1=PH⊤(HPH⊤+R)−1=HPH⊤+RPH⊤这就是卡尔曼增益的标准形式。(写成分数形式是为了方便理解)
观察分母,当测量噪声协方差 R 很大时,说明测量值不可靠,卡尔曼增益 K→0 ,更相信先验估计;当 R 很小时,说明测量值比较可靠,卡尔曼增益 K→H− ,更相信测量值。
一个完整的卡尔曼滤波器迭代流程#
1.预测阶段#
1.1 预测状态#
求先验估计(预测状态)x^k∣k−1,根据 k−1 时刻的后验估计(更新状态) x^k−1∣k−1 和系统的状态传播模型来预测 k 时刻的状态:
x^k∣k−1=Fkx^k−1∣k−1+Bkuk1.2 预测协方差矩阵#
求 Pk∣k−1,根据 k−1 时刻的后验协方差矩阵 Pk−1∣k−1 和系统的过程噪声协方差 Qk 来预测 k 时刻的协方差:
Pk∣k−1=FkPk−1∣k−1FkT+Qk2.更新阶段#
2.1 计算残差#
计算测量残差 yk,即当前测量值 zk 与估计测量值 z^k 之间的差异:
yk=zk−z^k=zk−Hkx^k∣k−12.2 计算卡尔曼增益#
计算卡尔曼增益 Kk,它决定了先验估计和测量值在更新中的权重:
SkKk=HkPk∣k−1Hk⊤+Rk=Pk∣k−1Hk⊤Sk−1=Pk∣k−1Hk⊤(HkPk∣k−1Hk⊤+Rk)−12.3 算后验估计#
更新之前的先验估计 x^k∣k−1,得到后验估计(更新后的状态)x^k∣k:
x^k∣k=x^k∣k−1+Kkyk2.4 更新协方差矩阵#
更新之前的先验协方差矩阵 Pk∣k−1,得到后验协方差矩阵(更新后的估计精度)Pk∣k:
Pk∣k=(I−KkHk)Pk∣k−13.迭代#
重复步骤1,2,直到达到所需的估计精度或处理完所有测量数据。
flowchart TD
A[预测阶段] --> B[预测状态]
B --> C[预测协方差矩阵]
C --> D[更新阶段]
D --> E[计算残差]
E --> F[计算卡尔曼增益 ]
F --> G[更新状态估计]
G --> H[更新协方差矩阵]
H --> A
恭喜你速通成功,快去看看卡尔曼滤波的代码实现吧!