卡尔曼滤波器
卡尔曼滤波器是贝叶斯滤波器在线性高斯系统下的一个特例。
假设线性高斯系统下运动模型和观测模型如下:
运动方程: x k = A k − 1 x k − 1 + v k + w k , k = 1 , ⋯ , K 观测方程: y k = C k x k + n k , k = 0 , ⋯ , K \begin{aligned}&\text{运动方程:}\boldsymbol{x}_{k}=\boldsymbol{A}_{k-1} \boldsymbol{x}_{k-1}+\boldsymbol {v}_{k}+\boldsymbol{w}_{k}, \quad k=1, \cdots, K \\&\text{观测方程:}\boldsymbol{y}_{k}=\boldsymbol{C}_{k} x_{k}+\boldsymbol{n}_{k}, \quad k=0, \cdots, K\end{aligned} 运动方程:xk=Ak−1xk−1+vk+wk,k=1,⋯,K观测方程:yk=Ckxk+nk,k=0,⋯,K
其中:
- x k ∈ R N \boldsymbol{x}_{k} \in \mathbb{R}^{N} xk∈RN 表示系统状态
- v k ∈ R N \boldsymbol{v}_{k} \in \mathbb{R}^{N} vk∈RN 表示输入
- w k ∈ R N \boldsymbol{w}_{k} \in \mathbb{R}^{N} wk∈RN 表示过程噪声, w k ∼ N ( 0 , Q k ) \boldsymbol{w}_{k} \sim N(\mathbf{0}, \boldsymbol{Q}_{k}) wk∼N(0,Qk)
- y k ∈ R M \boldsymbol{y}_{k} \in \mathbb{R}^{M} yk∈RM 表示测量
- n k ∈ R M \boldsymbol{n}_{k} \in \mathbb{R}^{M} nk∈RM 表示测量噪声, n k ∼ N ( 0 , R k ) \boldsymbol{n}_{k} \sim N(\mathbf{0}, \boldsymbol{R}_{k}) nk∼N(0,Rk)
- k k k 为时间下标,最大值为 K K K
用 ( ⋅ ) ˇ \check{(\cdot)} (⋅)ˇ 表示先验,用 ( ⋅ ) ^ \hat{(\cdot)} (⋅)^ 表示后验,卡尔曼滤波器则可以表示为:
预测
x ˇ k = A k − 1 x ^ k − 1 + v k \check{\boldsymbol{x}}_{k}=\boldsymbol{A}_{k-1} \hat{\boldsymbol{x}}_{k-1}+\boldsymbol{v}_{k} xˇk=Ak−1x^k−1+vk
P ˇ k = A k − 1 P ^ k − 1 A k − 1 T + Q k \check{\boldsymbol{P}}_{k}=\boldsymbol{A}_{k-1} \hat{\boldsymbol{P}}_{k-1} \boldsymbol{A}_{k-1}^{\mathrm{T}}+\boldsymbol{Q}_{k} Pˇk=Ak−1P^k−1Ak−1T+Qk
卡尔曼增益
K k = P ˇ k C k T ( C k P ˇ k C k T + R k ) − 1 \quad \boldsymbol{K}_{k}=\check{\boldsymbol{P}}_{k} \boldsymbol{C}_{k}^{\mathrm{T}}\left(\boldsymbol{C}_{k} \check{\boldsymbol{P}}_{k} \boldsymbol{C}_{k}^{\mathrm{T}}+\boldsymbol{R}_{k}\right)^{-1} Kk=PˇkCkT(CkPˇkCkT+Rk)−1
更新
x ^ k = x ˇ k + K k ( y k − C k x ˇ k ) ⏟ 更新量 \hat{\boldsymbol{x}}_{k}=\check{\boldsymbol{x}}_{k}+\boldsymbol{K}_{k} \underbrace{\left(\boldsymbol{y}_{k}-\boldsymbol{C}_{k} \check{\boldsymbol{x}}_{k}\right)}_{\text {更新量 }} x^k=xˇk+Kk更新量 (yk−Ckxˇk)
P ^ k = ( 1 − K k C k ) P ˇ k \hat{\boldsymbol{P}}_{k}=\left(\mathbf{1}-\boldsymbol{K}_{k} \boldsymbol{C}_{k}\right) \check{\boldsymbol{P}}_{k} P^k=(1−KkCk)Pˇk
通过贝叶斯推断推导
预测
由于是线性高斯系统,因此直接将 k − 1 k-1 k−1 时刻的后验分布通过线性运动模型传递,即可得到 k k k 时刻的先验:
x ˇ k = E [ x k ] = E [ A k − 1 x k − 1 + v k + w k ] = A k − 1 E [ x k − 1 ] ⏟ x k − 1 + v k + E [ w k ] ⏟ 0 = A k − 1 x ^ k − 1 + v k \begin{aligned}\check{\boldsymbol{x}}_{k} &=E\left[\boldsymbol{x}_{k}\right]\\&=E\left[\boldsymbol{A}_{k-1} \boldsymbol{x}_{k-1}+\boldsymbol{v}_{k}+\boldsymbol{w}_{k}\right] \\&=\boldsymbol{A}_{k-1} \underbrace{E\left[\boldsymbol{x}_{k-1}\right]}_{\boldsymbol{x}_{k-1}}+\boldsymbol{v}_{k}+\underbrace{E\left[\boldsymbol{w}_{k}\right]}_{0}\\&=\boldsymbol{A}_{k-1} \hat{\boldsymbol{x}}_{k-1}+\boldsymbol{v}_{k}\end{aligned} xˇk=E[xk]=E[Ak−1xk−1+vk+wk]=Ak−1xk−1 E[xk−1]+vk+0 E[wk]=Ak−1x^k−1+vk
P ˇ k = E [ ( x k − E [ x k ] ) ( x k − E [ x k ] ) T ] = E [ ( A k − 1 x k − 1 + v k + w k − A k − 1 x ^ k − 1 − v k ) × ( A k − 1 x k − 1 + v k + w k − A k − 1 x ^ k − 1 − v k ) T ] = A k − 1 E [ ( x k − 1 − x ^ k − 1 ) ( x k − 1 − x ^ k − 1 ) T ] ⏟ P ^ k − 1 A k − 1 T + E [ w k w k T ] ⏟ Q k = A k − 1 P ^ k − 1 A k − 1 T + Q k \begin{aligned}\check{\boldsymbol{P}}_{k}=& E\left[\left(\boldsymbol{x}_{k}-E\left[\boldsymbol{x}_{k}\right]\right)\left(\boldsymbol{x}_{k}-E\left[\boldsymbol{x}_{k}\right]\right)^{\mathrm{T}}\right] \\=& E\left[\left(\boldsymbol{A}_{k-1} \boldsymbol{x}_{k-1}+\boldsymbol{v}_{k}+\boldsymbol{w}_{k}-\boldsymbol{A}_{k-1} \hat{\boldsymbol{x}}_{k-1}-\boldsymbol{v}_{k}\right)\right.\\&\left.\times\left(\boldsymbol{A}_{k-1} \boldsymbol{x}_{k-1}+\boldsymbol{v}_{k}+\boldsymbol{w}_{k}-\boldsymbol{A}_{k-1} \hat{\boldsymbol{x}}_{k-1}-\boldsymbol{v}_{k}\right)^{\mathrm{T}}\right] \\=& \boldsymbol{A}_{k-1} \underbrace{E\left[\left(\boldsymbol{x}_{k-1}-\hat{\boldsymbol{x}}_{k-1}\right)\left(\boldsymbol{x}_{k-1}-\hat{\boldsymbol{x}}_{k-1}\right)^{\mathrm{T}}\right]}_{\hat{\boldsymbol{P}}_{k-1}} \boldsymbol{A}_{k-1}^{\mathrm{T}}+\underbrace{E\left[\boldsymbol{w}_{k} \boldsymbol{w}_{k}^{\mathrm{T}}\right]}_{\boldsymbol{Q}_{k}} \\=& \boldsymbol{A}_{k-1} \hat{\boldsymbol{P}}_{k-1} \boldsymbol{A}_{k-1}^{\mathrm{T}}+\boldsymbol{Q}_{k}\end{aligned} Pˇk====E[(xk−E[xk])(xk−E[xk])T]E[(Ak−1xk−1+vk+wk−Ak−1x^k−1−vk)×(Ak−1xk−1+vk+wk−Ak−1x^k−1−vk)T]Ak−1P^k−1 E[(xk−1−x^k−1)(xk−1−x^k−1)T]Ak−1T+Qk E[wkwkT]Ak−1P^k−1Ak−1T+Qk
更新
对于更新部分,我们将状态与 k k k 时刻的测量写成联合高斯分布的形式:
p ( x k , y k ∣ x ˇ 0 , v 1 : k , y 0 : k − 1 ) = N ( [ μ x μ y ] , [ Σ x x Σ x y Σ y x Σ y y ] ) = N ( [ x ˇ k C k x ˇ k ] , [ P ˇ k P ˇ k C k T C k P ˇ k C k P ˇ k C k T + R k ] ) \begin{aligned}p\left(\boldsymbol{x}_{k}, \boldsymbol{y}_{k} \mid \check{\boldsymbol{x}}_{0}, \boldsymbol{v}_{1: k}, \boldsymbol{y}_{0: k-1}\right) &=\mathcal{N}\left(\left[\begin{array}{c}\boldsymbol{\mu}_{x} \\\boldsymbol{\mu}_{y}\end{array}\right],\left[\begin{array}{cc}\boldsymbol{\Sigma}_{x x} & \boldsymbol{\Sigma}_{x y} \\\boldsymbol{\Sigma}_{y x} & \boldsymbol{\Sigma}_{y y}\end{array}\right]\right) \\&=\mathcal{N}\left(\left[\begin{array}{c}\check{\boldsymbol{x}}_{k} \\\boldsymbol{C}_{k} \check{\boldsymbol{x}}_{k}\end{array}\right],\left[\begin{array}{cc}\check{\boldsymbol{P}}_{k} & \check{\boldsymbol{P}}_{k} \boldsymbol{C}_{k}^{\mathrm{T}} \\\boldsymbol{C}_{k} \check{\boldsymbol{P}}_{k} & \boldsymbol{C}_{k} \check{\boldsymbol{P}}_{k} \boldsymbol{C}_{k}^{\mathrm{T}}+\boldsymbol{R}_{k}\end{array}\right]\right)\end{aligned} p(xk,yk∣xˇ0,v1:k,y0:k−1)=N([μxμy],[ΣxxΣyxΣxyΣyy])=N([xˇkCkxˇk],[PˇkCkPˇkPˇkCkTCkPˇkCkT+Rk])
由高斯推断可以直接得到 k k k 时刻的条件分布(即后验概率):
p ( x k ∣ x ˇ 0 , v 1 : k , y 0 : k ) = N ( μ x + Σ x y Σ y y − 1 ( y k − μ y ) ⏟ x ^ k , Σ x x − Σ x y Σ y y − 1 Σ y x ⏟ P ^ k ) p\left(\boldsymbol{x}_{k} \mid \check{\boldsymbol{x}}_{0}, \boldsymbol{v}_{1: k}, \boldsymbol{y}_{0: k}\right)=\mathcal{N}(\underbrace{\boldsymbol{\mu}_{x}+\boldsymbol{\Sigma}_{x y} \boldsymbol{\Sigma}_{y y}^{-1}\left(\boldsymbol{y}_{k}-\boldsymbol{\mu}_{y}\right)}_{\hat{\boldsymbol{x}}_{k}}, \underbrace{\boldsymbol{\Sigma}_{x x}-\boldsymbol{\Sigma}_{x y} \boldsymbol{\Sigma}_{y y}^{-1} \boldsymbol{\Sigma}_{y x}}_{\hat{\boldsymbol{P}}_{k}}) p(xk∣xˇ0,v1:k,y0:k)=N(x^k μx+ΣxyΣyy−1(yk−μy),P^k Σxx−ΣxyΣyy−1Σyx)
令 K k = Σ x y Σ y y − 1 = P ˇ k C k T ( C k P ˇ k C k T + R k ) − 1 \boldsymbol{K}_{k}=\boldsymbol{\Sigma}_{x y} \boldsymbol{\Sigma}_{y y}^{-1}=\check{\boldsymbol{P}}_{k} \boldsymbol{C}_{k}^{\mathrm{T}}\left(\boldsymbol{C}_{k} \check{\boldsymbol{P}}_{k} \boldsymbol{C}_{k}^{\mathrm{T}}+\boldsymbol{R}_{k}\right)^{-1} Kk=ΣxyΣyy−1=PˇkCkT(CkPˇkCkT+Rk)−1(卡尔曼增益),则有:
x ^ k = x ˇ k + K k ( y k − C k x ˇ k ) \hat{\boldsymbol{x}}_{k}=\check{\boldsymbol{x}}_{k}+\boldsymbol{K}_{k} \left(\boldsymbol{y}_{k}-\boldsymbol{C}_{k} \check{\boldsymbol{x}}_{k}\right) x^k=xˇk+Kk(yk−Ckxˇk)
P ^ k = ( 1 − K k C k ) P ˇ k \hat{\boldsymbol{P}}_{k}=\left(\mathbf{1}-\boldsymbol{K}_{k} \boldsymbol{C}_{k}\right) \check{\boldsymbol{P}}_{k} P^k=(1−KkCk)Pˇk
即证。