扩展卡尔曼滤波器

卡尔曼滤波器是最优的线性滤波器,在非线性系统中,没有办法用线性的状态空间方程表达,只能以如下形式表达:
Xk=f(Xk−1,Uk−1,Wk−1)Zk=h(Xk,Vk) \begin{align} X_k&=f(X_{k-1}, U_{k-1}, W_{k-1})\\ Z_k&=h(X_k, V_k) \end{align} Xk​Zk​​=f(Xk−1​,Uk−1​,Wk−1​)=h(Xk​,Vk​)​​
f、hf、hf、h为非线性函数,然而,正态分布的随机变量通过非线性系统后就不再是正态的了,故我们需要对其线性化,利用泰勒展开。
f(x)=f(x0)+∂f∂x(x−x0) f(x)=f(x_0)+\frac{\partial f}{\partial x}(x-x_0) f(x)=f(x0​)+∂x∂f​(x−x0​)
看得出来,若我们想要线性化一个系统,我们需要找到一个真实点,在其附近进行线性化。但系统有误差,我们没法知道真实点,则我们可以对于f(xk)f(x_k)f(xk​),在上一次的后验估计值x^k−1\hat x_{k-1}x^k−1​——即最佳估计值附近线性化:
Xk≈f(X^k−1,Uk−1)+Jk−1⏟∂f∂X∣X=X^k−1(Xk−1−X^k−1)+Wk−1 X_k \approx f(\hat X_{k-1}, U_{k-1}) + \underbrace {J_{k-1}}_{\frac{\partial f}{\partial X}|_{X=\hat X_{k-1}}}(X_{k-1} - \hat X_{k-1}) + W_{k-1} Xk​≈f(X^k−1​,Uk−1​)+∂X∂f​∣X=X^k−1​​Jk−1​​​(Xk−1​−X^k−1​)+Wk−1​
对于观测函数:
Zk≈h(X^k−)+Hk⏟∂h∂X∣X=X^k−(Xk−X^k−)+Vk Z_k \approx h(\hat X_k^-)+\underbrace{H_k}_{\frac{\partial h}{\partial X}|_{X=\hat X_k^-}}(X_k-\hat X_k^-) + V_k Zk​≈h(X^k−​)+∂X∂h​∣X=X^k−​​Hk​​​(Xk​−X^k−​)+Vk​
先验估计(使用非线性函数直接观测):
X^k−=f(X^k−1,Uk−1) \hat X_k^- = f(\hat X_{k-1}, U_k-1) X^k−​=f(X^k−1​,Uk​−1)
先验误差协方差:

首先计算误差:
ek−=Xk−X^k−①Xk≈f(X^k−1,Uk−1)+Jk−1(Xk−1−X^k−1)+Wk−1②X^k−=f(X^k−1,Uk−1)③⇒ek−≈Jk−1ek−1+Wk−1 \begin{align} e_k^- &= X_k - \hat X_k^- && ①\\ X_k &\approx f(\hat X_{k-1}, U_{k-1}) + J_{k-1}(X_{k-1} - \hat X_{k-1}) + W_{k-1} && ②\\ \hat X_k^- &= f(\hat X_{k-1}, U_{k-1}) && ③\\ \Rightarrow e_k^- &\approx J_{k-1}e_{k-1}+W_{k-1} \end{align} ek−​Xk​X^k−​⇒ek−​​=Xk​−X^k−​≈f(X^k−1​,Uk−1​)+Jk−1​(Xk−1​−X^k−1​)+Wk−1​=f(X^k−1​,Uk−1​)≈Jk−1​ek−1​+Wk−1​​​①②③​​
协方差传播:
Pk−=E[ek−ek−T]=E[(Jk−1ek−1+Wk−1)(Jk−1ek−1+Wk−1)T]=Jk−1E[ek−1ek−1T]Jk−1T+E[Wk−1Wk−1T]=Jk−1Pk−1Jk−1T+Qk−1 \begin{align} P_{k}^- &= E[e_k^-e_k^{-T}]\\ &=E[(J_{k-1}e_{k-1}+W_{k-1})(J_{k-1}e_{k-1}+W_{k-1})^T]\\ &=J_{k-1}E[e_{k-1}e_{k-1}^T]J_{k-1}^T+E[W_{k-1}W_{k-1}^T]\\ &=J_{k-1}P_{k-1}J_{k-1}^T+Q_{k-1} \end{align} Pk−​​=E[ek−​ek−T​]=E[(Jk−1​ek−1​+Wk−1​)(Jk−1​ek−1​+Wk−1​)T]=Jk−1​E[ek−1​ek−1T​]Jk−1T​+E[Wk−1​Wk−1T​]=Jk−1​Pk−1​Jk−1T​+Qk−1​​​
后验误差协方差:

后验估计误差:
ek=Xk−X^k=Xk−X^k−−Kk(Zk−h(X^k−))≈ek−−KkHkeK−−KkVk≈(I−KkHk)ek−−KkVk \begin{align} e_k &= X_k - \hat X_k\\ &= X_k - \hat X_k^--K_k(Z_k - h(\hat X_k^-))\\ &\approx e_k^- - K_kH_ke_K^--K_kV_k\\ &\approx(I-K_kH_k)e_k^--K_kV_k \end{align} ek​​=Xk​−X^k​=Xk​−X^k−​−Kk​(Zk​−h(X^k−​))≈ek−​−Kk​Hk​eK−​−Kk​Vk​≈(I−Kk​Hk​)ek−​−Kk​Vk​​​

则后验误差协方差:
Pk=E[ekekT]=E[((I−KkHk)ek−−KkVk)((I−KkHk)ek−−KkVk)T]=E[((I−KkHk)ek−−KkVk)(ek−T(I−HkTKkT)−KkTVkT)]=E[(I−KkHk)ek−ek−T(I−HkTKkT)−(I−KkHk)eK−KkTVkT−KkVkek−T(I−HkTKkT)+KkVkVkTKkT]=E[(I−KkHk)ek−ek−T(I−HkTKkT)]−E[(I−KkHk)eK−KkTVkT]−E[KkVkek−T(I−HkTKkT)]+E[KkVkVkTKkT]=(I−KkHk)Pk−(I−HkTKkT)+KkRkKkT⏟协方差更新的Joseph形式=Pk−−KkHkPk−−Pk−HkTKkT+Kk(HPk−HT+R)KkT \begin{align} P_k&=E[e_ke_k^T]\\ &=E[((I-K_kH_k)e_k^--K_kV_k)((I-K_kH_k)e_k^--K_kV_k)^T]\\ &=E[((I-K_kH_k)e_k^--K_kV_k)(e_k^{-T}(I-H_k^TK_k^T)-K_k^TV_k^T)]\\ &=E[(I-K_kH_k)e_k^-e_k^{-T}(I-H_k^TK_k^T)-(I-K_kH_k)e_K^-K_k^TV_k^T-K_kV_ke_k^{-T}(I-H_k^TK_k^T)+K_kV_kV_k^TK_k^T]\\ &=E[(I-K_kH_k)e_k^-e_k^{-T}(I-H_k^TK_k^T)]-E[(I-K_kH_k)e_K^-K_k^TV_k^T]-E[K_kV_ke_k^{-T}(I-H_k^TK_k^T)]+E[K_kV_kV_k^TK_k^T]\\ &=\underbrace{(I-K_kH_k)P_k^-(I-H_k^TK_k^T) + K_kR_kK_k^T}_{\text{协方差更新的Joseph形式}}\\ &=P_k^--K_kH_kP_k^--P_k^-H_k^TK_k^T+K_k(HP_k^-H^T+R)K_k^T \end{align} Pk​​=E[ek​ekT​]=E[((I−Kk​Hk​)ek−​−Kk​Vk​)((I−Kk​Hk​)ek−​−Kk​Vk​)T]=E[((I−Kk​Hk​)ek−​−Kk​Vk​)(ek−T​(I−HkT​KkT​)−KkT​VkT​)]=E[(I−Kk​Hk​)ek−​ek−T​(I−HkT​KkT​)−(I−Kk​Hk​)eK−​KkT​VkT​−Kk​Vk​ek−T​(I−HkT​KkT​)+Kk​Vk​VkT​KkT​]=E[(I−Kk​Hk​)ek−​ek−T​(I−HkT​KkT​)]−E[(I−Kk​Hk​)eK−​KkT​VkT​]−E[Kk​Vk​ek−T​(I−HkT​KkT​)]+E[Kk​Vk​VkT​KkT​]=协方差更新的Joseph形式(I−Kk​Hk​)Pk−​(I−HkT​KkT​)+Kk​Rk​KkT​​​=Pk−​−Kk​Hk​Pk−​−Pk−​HkT​KkT​+Kk​(HPk−​HT+R)KkT​​​
求解最优增益:

为了最小化PkP_kPk​的迹,对KkK_kKk​求偏导并令其为0:
∂tr(Pk)∂Kk=−2Pk−HkT+2Kk(HkPk−HT+R)=0 \begin{align} \frac{\partial tr(P_k)}{\partial K_k}=-2P_k^-H_k^T+2K_k(H_kP_k^-H^T+R)=0 \end{align} ∂Kk​∂tr(Pk​)​=−2Pk−​HkT​+2Kk​(Hk​Pk−​HT+R)=0​​
则最优卡尔曼增益为:
Kk=Pk−HkT(HPk−HT+R)−1 K_k=P_k^-H_k^T(HP_k^-HT+R)^{-1} Kk​=Pk−​HkT​(HPk−​HT+R)−1
状态更新:

利用最优增益更新状态:
X^k=X^k−+Kk[Zk−h(X^h−)] \hat X_k = \hat X_k^- + K_k[Z_k-h(\hat X_h^-)] X^k​=X^k−​+Kk​[Zk​−h(X^h−​)]
若KkK_kKk​较大,更信任测量值,测量噪声小

若KkK_kKk​较小,更信任预测值,预测误差小

协方差更新:

将最优增益带入PkP_kPk​,可得:
Pk=(I−KkHk)Pk− \begin{align} P_k &= (I-K_kH_k)P_k^- \end{align} Pk​​=(I−Kk​Hk​)Pk−​​​

完整流程

初始化:

X^0=E[X0]P0=Cov(X0) \begin{align} \hat X_0 &= E[X_0]\\ P_0 &= Cov(X_0) \end{align} X^0​P0​​=E[X0​]=Cov(X0​)​​

预测阶段:

状态雅可比矩阵:
Jk−1=∂f∂X∣X=X^k−1 J_{k-1} = \frac{\partial f}{\partial X}|_{X=\hat X_{k-1}} Jk−1​=∂X∂f​∣X=X^k−1​​
状态预测(先验估计):
X^k−=f(X^k−1+Uk−1) \hat X_k^-=f(\hat X_{k-1}+U_{k-1}) X^k−​=f(X^k−1​+Uk−1​)

误差协方差预测(先验误差协方差):
Pk−=Jk−1Pk−1Jk−1T+Qk−1 P_k^-=J_{k-1}P_{k-1}J_{k-1}^T + Q_{k-1} Pk−​=Jk−1​Pk−1​Jk−1T​+Qk−1​

更新阶段:

观测雅可比矩阵:
Hk=∂h∂X∣X=X^k− H_k = \frac {\partial h}{\partial X}|_{X=\hat X_k^-} Hk​=∂X∂h​∣X=X^k−​​

卡尔曼增益:
​ Kk=Pk−HTHPk−HT+RK_k=\frac{P_k^-H^T}{HP_k^-H^T+R}Kk​=HPk−​HT+RPk−​HT​

状态更新(后验估计):
​ X^k=X^k−+Kk(Zk−h(X^k−))\hat X_k = \hat X_k^-+K_k(Z_k-h(\hat X_k^-))X^k​=X^k−​+Kk​(Zk​−h(X^k−​))

误差协方差更新:
​ Pk=(I−KkH)Pk−P_k = (I-K_kH)P_k^-Pk​=(I−Kk​H)Pk−​
或Joseph\text{Joseph}Joseph形式:
Pk=(I−KkH)Pk−(I−KkH)T+KkRKkT. P_k = (I - K_k H)P_k^-(I - K_k H)^T + K_k R K_k^T. Pk​=(I−Kk​H)Pk−​(I−Kk​H)T+Kk​RKkT​.

Logo

北京人形旗下天工造物具身智能开源社区,聚焦具身天工与慧思开物两大平台

更多推荐