机器学习-线性动态系统-卡尔曼滤波(Kalman Filtering)
Dynamic Model
观测值:设有一系列观测值 y 1 , y 2 ⋯ y n y_1,y_2 \cdots y_n y1,y2⋯yn,不可以互换顺序,观测之间是相关的,有联系的。
状态空间模型(state-based-model):引入隐状态 x 1 , x 2 ⋯ x n x_1,x_2 \cdots x_n x1,x2⋯xn.在给定系列隐状态的条件下,观测之间变为相互独立
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gCnjpbgV-1636796568836)(C:\Users\nth12\AppData\Roaming\Typora\typora-user-images\image-20211112171837873.png)]
- 绿色的箭头指示了 p ( x t ∣ x t − 1 ) − − t r a n s s i s o n p r o b p(x_t | x_{t-1})--transsison prob p(xt∣xt−1)−−transsisonprob
- 红色箭头指示了 p ( y t ∣ x t ) − − m e a s u r e m e n t p r o b , e m i s s i o n p r o b p(y_t|x_t)--measurement prob,emission prob p(yt∣xt)−−measurementprob,emissionprob
- 描述HMM: p ( x 1 ) − − i n i t p r o b , p ( x t ∣ x t − 1 ) , p ( y t ∣ x t ) ∀ t p(x_1)--init \ prob,p(x_t|x_{t-1}),p(y_t|x_t)\quad \forall t p(x1)−−init prob,p(xt∣xt−1),p(yt∣xt)∀t
离散状态模型(Discrete state Dynamic Model)
-
HMM(隐马尔可夫模型)-所有隐状态即 x t x_t xt离散
-
p ( x t ∣ x t − 1 ) p(x_t | x_{t-1}) p(xt∣xt−1)由A矩阵给出,其中A矩阵满足
A x t − 1 , x t = p ( x t ∣ x t − 1 ) A_{x_{t-1},x_t} = p(x_t|x_{t-1}) Axt−1,xt=p(xt∣xt−1) -
p ( y t ∣ x t ) p(y_t|x_t) p(yt∣xt)—任何形式。 y t y_t yt可以是离散的,可以是连续的。
-
p ( x 1 ) p(x_1) p(x1)由 π \pi π向量确定。
线性+高斯噪音模型(Linear Guass Dynamic Model)
-
卡尔曼滤波
-
p ( x t ∣ x t − 1 ) p(x_t|x_{t-1}) p(xt∣xt−1)转移概率不能用矩阵表示,使用函数表示
p ( x t ∣ x t − 1 ) ∼ N ( A x t − 1 + B , Q ) p(x_t|x_{t-1}) \sim N(Ax_{t-1}+B,Q) p(xt∣xt−1)∼N(Axt−1+B,Q) -
p ( y t ∣ x t ) p(y_t|x_{t}) p(yt∣xt)也是线性的,使用函数表示:
p ( y t ∣ x t ) ∼ N ( H x t − 1 + C , R ) p(y_t|x_t) \sim N(Hx_{t-1}+C,R) p(yt∣xt)∼N(Hxt−1+C,R) -
p ( x 1 ) ∼ N ( μ 0 , σ 0 ) p(x_1) \sim N(\mu_0,\sigma_0) p(x1)∼N(μ0,σ0)
非线性,非高斯的动态模型
-
粒子滤波
-
p ( x t ∣ x t − 1 ) p(x_t|x_{t-1}) p(xt∣xt−1)任意函数 f ( x t − 1 ) f(x_{t-1}) f(xt−1)
-
p ( y t ∣ x t ) ∼ N ( H x t − 1 + C , R ) p(y_t|x_t) \sim N(Hx_{t-1}+C,R) p(yt∣xt)∼N(Hxt−1+C,R)也是任意函数
-
p ( x 1 ) p(x_1) p(x1)也是任意函数
基本任务
- eval:求解 p ( y 1 , y 2 , ⋯ y n ) p(y_1,y_2,\cdots y_n) p(y1,y2,⋯yn)
- parameter learning: a r g m a x θ l o g p ( y 1 , ⋯ y n ∣ θ ) argmax_\theta \ logp(y_1,\cdots y_n | \theta) argmaxθ logp(y1,⋯yn∣θ)
- state decoding: a r g m a x x 1 , x 2 ⋯ x n p ( x 1 , x 2 ⋯ x n ∣ y 1 , y 2 , ⋯ y n ) argmax_{x_1,x_2\cdots x_n} \ p(x_1,x_2 \cdots x_n | y_1,y_2,\cdots y_n) argmaxx1,x2⋯xn p(x1,x2⋯xn∣y1,y2,⋯yn)
- filtering:求解 p ( x t ∣ y 1 ⋯ y t ) p(x_t | y_1 \cdots y_t) p(xt∣y1⋯yt)
卡尔曼滤波
动机
- 飞机、雷达的实时定位和参数测量都受到各种各样的随机干扰,要想正确的得到实时状态参数,只能根据观测到的信号
来预测真实状态 - 最优估计问题:希望估值误差尽可能小,产生了最优估计问题
- 目标:无偏性,估计方差最小,实时性
概述
-
转移概率 p ( x t ∣ x t − 1 ) = N ( A x t − 1 + B , Q ) p(x_t|x_{t-1}) = N(Ax_{t-1}+B,Q) p(xt∣xt−1)=N(Axt−1+B,Q)
- 这等价于: p ( x t ∣ x t − 1 ) = A x t − 1 + B + w p(x_t|x_{t-1}) = Ax_{t-1} + B + w p(xt∣xt−1)=Axt−1+B+w,其中 w ∼ N ( 0 , Q ) w \sim N(0,Q) w∼N(0,Q)
-
measurement prob: p ( y t ∣ x t ) = N ( H x t + C , R ) p(y_t|x_t) = N(Hx_t+C,R) p(yt∣xt)=N(Hxt+C,R)
- 这等价于: p ( y t ∣ x t ) = H x t + C + V , V ∼ N ( 0 , R ) p(y_t | x_t) = Hx_t + C + V,V \sim N(0,R) p(yt∣xt)=Hxt+C+V,V∼N(0,R)
-
所以其中的参数有:A,B,Q,R,C,H;
-
现在假设 A = I A=I A=I, B = ( 2 , 2 ) T B = (2,2)^T B=(2,2)T
- p ( y t ∣ x t ) = I x t + V , V ∼ N ( 0 , I ) p(y_t|x_t) = Ix_t+ V,V \sim N(0,I) p(yt∣xt)=Ixt+V,V∼N(0,I)
- p ( x t ∣ x t − 1 ) = I x t − 1 + ( 2 , 2 ) T + w p(x_t | x_{t-1}) = Ix_{t-1} + (2,2)^T + w p(xt∣xt−1)=Ixt−1+(2,2)T+w, w ∼ N ( 0 , I ) w \sim N(0,I) w∼N(0,I)
- 假设 x 1 = ( 0 , 0 ) T x_1 = (0,0)^T x1=(0,0)T
- 则
x
2
x_2
x2的位置为以(2,2)为中心的,I为方差的范围内。
- 其中 x 2 x_2 x2位于 ( 2 , 2 ) (2,2) (2,2)的概率密度是最大的,因为 w = 0 w = 0 w=0的对应概率密度最大。
- 再由式子 p ( y 2 ∣ x 2 ) p(y_2 | x_2) p(y2∣x2)确定 y 2 y_2 y2的取值, y 2 y_2 y2取值范围为 x 2 x_2 x2为中心,方差为I的高斯分布。
- 现在欲求解的问题是,我们观测到了 y 1 , y 2 ⋯ y t y_1,y_2\cdots y_t y1,y2⋯yt,欲求 p ( x t ∣ y 1 ⋯ y t ) p(x_t | y_1 \cdots y_t) p(xt∣y1⋯yt)
-
物体运动的实例:设一小车在直线上运动,其加速度 a ∼ N ( 0 , σ ) a \sim N(0,\sigma) a∼N(0,σ)
- 这里的 y t y_t yt在t时刻的位置(数轴上的位置坐标)
- 状态向量 x ‾ t = ( x t , x t , ) T \overline x_t = (x_t,x_t^,)^T xt=(xt,xt,)T,(与前面的向量进行区分)其中前一个分量表示汽车的真实位置,后一个分量表示汽车的真实速度
- 则有: x t = x t − 1 + x t − 1 , Δ t + 1 2 a Δ t 2 x_t = x_{t-1} + x_{t-1}^,\Delta t+\frac{1}{2}a\Delta t^2 xt=xt−1+xt−1,Δt+21aΔt2, x t , = x t − 1 , + a Δ t x^,_{t} = x^,_{t-1} + a \Delta t xt,=xt−1,+aΔt.使用这两个等式来表示卡尔曼滤波中的参数。
基本假设
方程
x
k
=
A
x
k
−
1
+
B
u
k
−
1
+
q
k
−
1
y
k
=
H
x
k
+
r
k
x_k = Ax_{k-1}+Bu_{k - 1}+q_{k-1} \\ y_k = Hx_k + r_k
xk=Axk−1+Buk−1+qk−1yk=Hxk+rk
-
其中 x k x_k xk是k时刻系统的状态; y k y_k yk是k时刻系统的测量值。 u k u_k uk是k时刻对系统的控制量
-
A,B是系统参数,对于多模型系统,他们为矩阵;
-
H为测量系统的参数,对于多测量系统,H为矩阵
-
q k q_k qk和 r k r_k rk表示过程和测量的噪声。
-
被假设称高斯白噪声( w h i t e g u a s s n o i s e white\ guass \ noise white guass noise),协方差矩阵分别为 Q , R Q,R Q,R
-
假设他们不随系统状态的变化而变化
-
高斯白噪声:功率谱密度服从均匀分布,幅度服从高斯分布
-
若N(t)为一个具有零均值的平稳随机过程,其功率谱密度均匀分布在 − ∞ , + ∞ -\infty,+\infty −∞,+∞ 的整个频率区间,即
S N ( ω ) = 1 2 N 0 S_N(\omega) = \frac{1}{2}N_0 SN(ω)=21N0
其中 N 0 N_0 N0为一正实数。
-
- 以最小均方误差为估计的最佳准则
- 利用前一时刻的估计值和现时刻的观测值来更新对状态变量的估计
原理推导
设线性系统的状态方程为
x
k
=
A
x
k
−
1
+
B
u
k
−
1
+
q
k
−
1
x_k = Ax_{k-1} + Bu_{k - 1}+q_{k-1}
xk=Axk−1+Buk−1+qk−1
- 其中 x x x为系统的状态变量,大小为 ( n , 1 ) (n,1) (n,1)。
- A为转换矩阵,维度为(n,n)。
- u为系统输入,维度为(k,1).
- B是将输入转换为状态的矩阵,维度为(n,k)。
- q为随机变量,系统噪声。假设服从分布 q ∼ N ( 0 , Q ) q\sim N(0,Q) q∼N(0,Q)
在概率论中,前面预测的结果称为先验,后面测量出的结果称为后验。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LqqOrZWW-1636796568843)(C:\Users\nth12\AppData\Roaming\Typora\typora-user-images\image-20211112162753126.png)]
测量值为后验,给出的公式为
z
k
=
H
x
k
+
r
k
z_k = Hx_k +r_k
zk=Hxk+rk
其中
z
k
z_k
zk的维度为(m,1),矩阵H的维度为(m,n),随机变量r为测量噪声,
r
∼
N
(
0
,
R
)
r \sim N(0,R)
r∼N(0,R),且与Q互相独立。
存中…(img-LqqOrZWW-1636796568843)]
测量值为后验,给出的公式为
z
k
=
H
x
k
+
r
k
z_k = Hx_k +r_k
zk=Hxk+rk
其中
z
k
z_k
zk的维度为(m,1),矩阵H的维度为(m,n),随机变量r为测量噪声,
r
∼
N
(
0
,
R
)
r \sim N(0,R)
r∼N(0,R),且与Q互相独立。