Kalman filter算法介绍及Python实现
Kalman filter算法介绍及Python实现一、算法思路1.1 Kalman filter简介1.2 算法推导二、Python复现三、参考文章一、算法思路1.1 Kalman filter简介卡尔曼滤波(Kalman filtering)是一种利用线性系统状态方程,通过系统输入输出观测数据,对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响,所以最优估计也可看作是滤波
Kalman filter算法介绍及Python实现
一、算法思路
1.1 Kalman filter简介
卡尔曼滤波(Kalman filtering)是一种利用线性系统状态方程,通过系统输入输出观测数据,对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响,所以最优估计也可看作是滤波过程。
数据滤波是去除噪声还原真实数据的一种数据处理技术,Kalman滤波在测量方差已知的情况下能够从一系列存在测量噪声的数据中,估计动态系统的状态。由于它便于计算机编程实现,并能够对现场采集的数据进行实时的更新和处理,Kalman滤波是目前应用最为广泛的滤波方法,在通信,导航,制导与控制等多领域得到了较好的应用。
1.2 算法推导
由于卡尔曼滤波是建立在三种状态下的预测模型,即现在我们现定义三种状态:
- 现实状态(1),也就是物体真实的轨迹
- 测量状态(2),也就是传感器返回的数值
- 预测状态(3),根据能动方程推论状态
{ X t = F t X t + B t U t + ω t ( 1 ) Z t = H t X t + V t ( 2 ) X ^ t = F t X ^ t − 1 + B t + μ t ( 3 ) \begin{cases} X_t = F_tX_t+B_tU_t+\omega_t \quad\quad\quad(1) \\Z_t = H_tX_t+V_t\quad\quad\quad\quad\quad\ \ \ \ \ (2) \\\hat {X}_{t} = F_t\hat {X}_{t-1}+B_t+\mu_t \quad\quad\quad(3) \end{cases} ⎩⎪⎨⎪⎧Xt=FtXt+BtUt+ωt(1)Zt=HtXt+Vt (2)X^t=FtX^t−1+Bt+μt(3)
其中(1)(2),大概就是位移随着时间的变化具有一定的线性关系, f ( x ) , h ( x ) f(x),h(x) f(x),h(x)就是这种抽象关系,可以故且看作一个基于运动方程的线性表出:
{ X t = f ( x t − 1 , μ t ) + ω t Z t , j = h ( y j , x t ) + v t , j \begin{cases} X_t = f(x_{t-1},\mu_t)+\omega_t \\Z_{t,j} = h(y_j,x_t)+v_{t,j} \end{cases} {Xt=f(xt−1,μt)+ωtZt,j=h(yj,xt)+vt,j
- x t , y j x_t,y_j xt,yj为 t t t时刻的状态向量,包括了相机位姿、路标坐标等信息,也可能有速度、朝向等信息(方便后续推断,将 y j y_j yj路标并入到状态里);
- μ t \mu_t μt 为运动测量值,如加速度,转向等等;
- F t F_t Ft 为状态转换方程,将 t − 1 t-1 t−1时刻的状态转换至 t t t时刻的状态;
- B t B_t Bt 是控制输入矩阵,将运动测量值 μ t \mu_t μt的作用映射到状态向量上;
- ω t \omega_t ωt是预测的高斯噪声,其均值为0,协方差矩阵为 Q t Q_t Qt。
其中(3),用一个在解释卡尔曼滤波时最常用的一维例子『小车追踪』,如下图所示:
由初中物理可得:
{ x = x t − 1 + v t × Δ t + 1 2 a Δ t 2 v t = v t − 1 + a t \begin{cases} x = x_{t-1}+v_t\times \Delta t+\frac{1}{2}a\Delta t^2 \\v_t = v_{t-1}+at \end{cases} {x=xt−1+vt×Δt+21aΔt2vt=vt−1+at
写成向量的形式就是:
[ x t v t ] = [ 1 Δ t 0 1 ] ∗ [ x t − 1 v t − 1 ] + [ 1 2 Δ t 2 Δ t ] ∗ a \left[ \begin{matrix} x_t\\v_t \end{matrix} \right] = \left[ \begin{matrix} 1& \Delta t \\ 0&1 \end{matrix} \right]* \left[ \begin{matrix} x_{t-1}\\v_{t-1} \end{matrix} \right]+ \left[ \begin{matrix} \frac{1}{2}\Delta t^2\\\Delta t \end{matrix} \right]*a [xtvt]=[10Δt1]∗[xt−1vt−1]+[21Δt2Δt]∗a
这时候你把公式等价化:
a ∼ μ t a \sim \mu_t a∼μt [ 1 Δ t 0 1 ] & & [ 1 2 Δ t 2 Δ t ] ∼ F t , B t \left[ \begin{matrix} 1& \Delta t \\ 0&1 \end{matrix} \right]\&\&\left[ \begin{matrix} \frac{1}{2}\Delta t^2\\\Delta t \end{matrix} \right]\sim F_t,B_t [10Δt1]&&[21Δt2Δt]∼Ft,Bt
就会发现动力学的预测公式就是这么个东西 X ^ t ∣ t − 1 = F t X ^ t − 1 + B t + μ t \hat {X}_{t|t-1} = F_t\hat {X}_{t-1}+B_t+\mu_t \quad\quad\quad X^t∣t−1=FtX^t−1+Bt+μt
- F t F_t Ft 为状态转换方程,将 t − 1 t-1 t−1时刻的状态转换至 t t t时刻的状态;
- B t B_t Bt 是控制输入矩阵,将运动测量值 μ t \mu_t μt的作用映射到状态向量上;
- X ^ t ∣ t − 1 \hat {X}_{t|t-1} X^t∣t−1 代表由 X t − 1 X_{t-1} Xt−1预测所来的 X t X_t Xt的值,为了和下面传感器的返回数值 X ^ t \hat X_t X^t区分就用这个来代替
用『式子(1)-(3)』:
X t − X ^ t ∣ t − 1 = F t ( X t − 1 − X ^ t − 1 ) + ω t X_t-\hat{X}_{t|t-1} = F_t(X_{t-1}-\hat X_{t-1})+\omega_t Xt−X^t∣t−1=Ft(Xt−1−X^t−1)+ωt
于是: P t ∣ t − 1 = E [ ( X t − X ^ t ∣ t − 1 ) ( X t − X ^ t ∣ t − 1 ) T ] … = F t P t − 1 F t T + Q t P_{t|t-1} = E[(X_t-\hat X_{t|t-1})(X_t-\hat X_{t|t-1})^T]\\\dots\\=F_tP_{t-1}F_t^T+Q_t Pt∣t−1=E[(Xt−X^t∣t−1)(Xt−X^t∣t−1)T]…=FtPt−1FtT+Qt
可以看到,经过预测更新,协方差矩阵P变大了。这是因为状态转换并不完美,而且运动测量值含有噪声,具有较大的不确定性。
预测更新实际上相当于“加法”:将当前状态转换到下一时刻(并增加一定不确定性),再把外界的干扰(运动测量值)叠加上去(又增加了一点不确定性)。
得到了测量值 X ^ \hat X X^,那么我们就可以对状态向量进行测量更新了,对应的公式为:
{ K t = P t ∣ t − 1 ∗ H t T H t P t ∣ t − 1 H T + R t P t = P t ∣ t − 1 − K t H t P t ∣ t − 1 \begin{cases} K_t = \frac{P_{t|t-1}*H_t^T}{H_tP_{t|t-1}H^T+R_t} \\P_t=P_{t|t-1}-K_tH_tP_{t|t-1} \end{cases} {Kt=HtPt∣t−1HT+RtPt∣t−1∗HtTPt=Pt∣t−1−KtHtPt∣t−1
可得:
X ^ t = X ^ t ∣ t − 1 + K t ( Z t − H t X ^ t ∣ t − 1 ) \hat X_t = \hat X_{t|t-1}+K_t(Z_t-H_t \hat X_{t|t-1}) X^t=X^t∣t−1+Kt(Zt−HtX^t∣t−1)
这一串的公式推导真的很麻烦,具体可以看知乎大佬的『卡尔曼滤波:从入门到精通』
二、Python复现
我这里是用的Github上的一个开源项目,具体的应用和实例文本Notebook链接:Asll’s Heywhale
三、参考文章
[1] 知乎:卡尔曼滤波:从入门到精通
[2] Github:pykalman
[3] CSDN:关于卡尔曼滤波和卡尔曼平滑关系的理解
[4] CSDN:一文图解卡尔曼滤波(Kalman Filter)
更多推荐
所有评论(0)