卡尔曼滤波
概述
卡尔曼滤波(Kalman filtering,KF)是一种利用线性系统状态方程,通过系统输入输出观测数据,对系统状态进行最优估计的算法。由于建模的不确定性导致系统本身存在过程噪声,以及利用传感器进行观测存在测量噪声,假设两者均服从正态分布,卡尔曼滤波在线性状态空间表示的基础上对有噪声的输入和观测信号进行处理,求取系统状态最优值。所以卡尔曼滤波适用于线性、离散和有限维系统。
卡尔曼滤波是目前应用最为广泛的滤波方法,在通信,导航,制导与控制等多领域得到了较好的应用。
公式及含义
假设线性系统的状态方程为:
x
k
=
A
x
k
−
1
+
B
u
k
−
1
+
w
k
−
1
x_k = Ax_{k-1}+Bu_{k-1}+w_{k-1}
xk=Axk−1+Buk−1+wk−1
其中,
- x k x_k xk和 x k − 1 x_{k-1} xk−1分别是k时刻和k-1时刻的状态值;
- A A A是状态转移矩阵;
- B B B是控制输入矩阵;
- u k − 1 u_{k-1} uk−1是k-1时刻的输入;
- w k − 1 w_{k-1} wk−1是k-1时刻的过程噪声,服从期望为 0 0 0,协方差矩阵为 Q Q Q 的正态分布(高斯分布)。
观测方程为:
z
k
=
H
x
k
+
v
k
z_k = Hx_k+v_k
zk=Hxk+vk
其中,
- x k x_k xk分别是k时刻的状态值;
- z k z_k zk是k时刻的测量值(观测值)
- v k v_k vk是k时刻的测量噪声(观测噪声),服从期望为 0 0 0,协方差矩阵为 R R R 的正态分布(高斯分布);
- H是状态观测矩阵。
卡尔曼滤波器时间更新方程:
x
^
k
−
=
A
x
^
k
−
1
+
B
u
k
−
1
P
k
−
=
A
P
k
−
1
A
T
+
Q
\hat{x}^{-}_{k} = A\hat{x}_{k-1} + Bu_{k-1} \\ P^{-}_{k} = AP_{k-1}A^{T} + Q
x^k−=Ax^k−1+Buk−1Pk−=APk−1AT+Q
卡尔曼滤波器状态更新方程:
K
k
=
P
k
−
H
T
H
P
k
−
H
T
+
R
x
^
k
=
x
^
k
−
+
K
k
(
z
k
−
H
x
^
k
−
)
P
k
=
(
I
−
K
k
H
)
P
k
−
K_k = \frac{P^{-}_{k}H^T}{HP^{-}_{k}H^{T}+R} \\ \hat{x}_{k} = \hat{x}^{-}_{k} + K_{k}(z_{k}-H\hat{x}^{-}_{k}) \\ P_{k} = (I-K_{k}H)P^{-}_{k}
Kk=HPk−HT+RPk−HTx^k=x^k−+Kk(zk−Hx^k−)Pk=(I−KkH)Pk−
其中,
- x ^ k \hat{x}_k x^k和 x ^ k − 1 \hat{x}_{k-1} x^k−1分别是k时刻和k-1时刻的后验估计值,即更卡尔曼滤波后的最优估计值;
- x ^ k − \hat{x}^{-}_k x^k−是k时刻的先验估计值,是卡尔曼滤波的中间结果之一;
- P k P_{k} Pk和 P k − 1 P_{k-1} Pk−1是k时刻和k-1时刻的后验估计协方差(即 x ^ k \hat{x}_k x^k和 x ^ k − 1 \hat{x}_{k-1} x^k−1的协方差),表示后验估计的不确定度;
- P k − P^{-}_{k} Pk−是k时刻的先验估计协方差,是卡尔曼滤波的中间结果之一;
- K k K_k Kk是卡尔曼增益系数,范围为;
- ( z k − H x ^ k − ) (z_{k}-H\hat{x}^{-}_{k}) (zk−Hx^k−)表示实际观测和预测观测的残差,和卡尔曼增益一起修正先验估计,得到后验估计值。
扩展卡尔曼滤波
概述
扩展卡尔曼滤波(Eetend Kalman Filtering, EKF)是对非线性系统状态方程进行线性化,再应用经典卡尔曼滤波算法对系统状态进行最优估计。
公式及含义
假非线性系统的状态方程为:
x
k
=
f
(
x
k
−
1
,
u
k
−
1
,
w
k
−
1
)
x_k = f(x_{k-1},u_{k-1},w_{k-1})
xk=f(xk−1,uk−1,wk−1)
其中,
- x k x_k xk和 x k − 1 x_{k-1} xk−1分别是k时刻和k-1时刻的状态值;
- u k − 1 u_{k-1} uk−1是k-1时刻的输入;
- w k − 1 w_{k-1} wk−1是k-1时刻的过程噪声,服从期望为 0 0 0,协方差矩阵为 Q Q Q 的正态分布(高斯分布)。
观测方程为:
z
k
=
h
(
x
k
,
v
k
)
z_k = h(x_k,v_k)
zk=h(xk,vk)
其中,
- x k x_k xk分别是k时刻的状态值;
- z k z_k zk是k时刻的测量值(观测值)
- v k v_k vk是k时刻的测量噪声(观测噪声),服从期望为 0 0 0,协方差矩阵为 R R R 的正态分布(高斯分布);
由于卡尔曼滤波只能应用于线性系统,需要对非线性系统进行线性化。
在
x
=
x
^
k
−
1
x=\hat{x}_{k-1}
x=x^k−1处进行线性化后的状态方程为:
x
k
=
x
~
k
+
A
(
x
k
−
1
−
x
^
k
−
1
)
+
W
k
w
k
−
1
x_k = \widetilde{x}_{k}+A(x_{k-1}-\hat{x}_{k-1})+W_kw_{k-1}
xk=x
k+A(xk−1−x^k−1)+Wkwk−1
其中,
- x ~ k = f ( x ^ k − 1 , u k − 1 , w k − 1 ) \widetilde{x}_{k}=f(\hat{x}_{k-1},u_{k-1},w_{k-1}) x k=f(x^k−1,uk−1,wk−1),由于 w k − 1 w_{k-1} wk−1未知,一般简化设置为0,所以 x ~ k = f ( x ^ k − 1 , u k − 1 , 0 ) \widetilde{x}_{k}=f(\hat{x}_{k-1},u_{k-1},0) x k=f(x^k−1,uk−1,0)
- A = ∂ f ∂ x ∣ x ^ k − 1 , u k − 1 A = \frac{\partial f}{\partial x}|_{\hat{x}_{k-1},u_{k-1}} A=∂x∂f∣x^k−1,uk−1
- W k = ∂ f ∂ w ∣ x ^ k − 1 , u k − 1 W_k = \frac{\partial f}{\partial w}|_{\hat{x}_{k-1},u_{k-1}} Wk=∂w∂f∣x^k−1,uk−1
在
x
=
x
~
k
x=\widetilde{x}_{k}
x=x
k处进行线性化后的观测方程为:
z
k
=
z
~
k
+
H
k
(
x
k
−
x
~
k
)
+
V
k
v
k
z_k = \widetilde{z}_k+H_k(x_k-\widetilde{x}_{k})+V_kv_{k}
zk=z
k+Hk(xk−x
k)+Vkvk
其中,
- z ~ k = h ( x ~ k , v k ) \widetilde{z}_{k}=h(\widetilde{x}_{k},v_{k}) z k=h(x k,vk),由于 v k v_{k} vk未知,一般简化设置为0,所以 z ~ k = h ( x ~ k , 0 ) \widetilde{z}_{k}=h(\widetilde{x}_{k},0) z k=h(x k,0)
- H k = ∂ h ∂ x ∣ x ~ k H_k = \frac{\partial h}{\partial x}|_{\widetilde{x}_{k}} Hk=∂x∂h∣x k
- V k = ∂ h ∂ v ∣ x ~ k V_k = \frac{\partial h}{\partial v}|_{\widetilde{x}_{k}} Vk=∂v∂h∣x k
扩展卡尔曼滤波器时间更新方程:
x
^
k
−
=
f
(
x
^
k
−
1
,
u
k
−
1
,
0
)
P
k
−
=
A
P
k
−
1
A
T
+
W
k
Q
k
−
1
W
k
T
\hat{x}^{-}_{k} = f(\hat{x}_{k-1},u_{k-1},0) \\ P^{-}_{k} = AP_{k-1}A^{T} + W_{k}Q_{k-1}W_{k}^T
x^k−=f(x^k−1,uk−1,0)Pk−=APk−1AT+WkQk−1WkT
扩展卡尔曼滤波器状态更新方程:
K
k
=
P
k
−
H
k
T
H
k
P
k
−
H
k
T
+
V
k
R
k
V
k
T
x
^
k
=
x
^
k
−
+
K
k
(
z
k
−
h
(
x
^
k
−
,
0
)
)
P
k
=
(
I
−
K
k
H
k
)
P
k
−
K_k = \frac{P^{-}_{k}H_k^T}{H_kP^{-}_{k}H_k^{T}+V_{k}R_{k}V_{k}^T} \\ \hat{x}_{k} = \hat{x}^{-}_{k} + K_{k}(z_{k}-h(\hat{x}^{-}_{k},0)) \\ P_{k} = (I-K_{k}H_k)P^{-}_{k}
Kk=HkPk−HkT+VkRkVkTPk−HkTx^k=x^k−+Kk(zk−h(x^k−,0))Pk=(I−KkHk)Pk−
Python代码实现
参考pykalman官方网站:https://pykalman.github.io/