一谈到贝叶斯滤波,就开始联系到各种随机过程、概率密度函数等等,那些曾经上课都听不进去的东西,这里能讲清楚吗?我自己也会有这个疑惑,不过要看懂贝叶斯滤波原理还是需要一定基础的。这篇博客,我会结合自己的理解尽量讲得通俗、方便理解一点。
这篇博客大部分参考了b站视频:忠厚老实的王大头
感谢up主做的视频
一、先验知识
1. 随机过程与概率论
两门都是大学学过的课程,那它们之间到底有什么关系呢?呃…其实大学期间自己也没理解,笑cry。
- 随机过程:研究的随机变量之间不独立,由于不独立,无法做随机实验,难度较高。
- 概率论:研究的随机变量之间需要独立,比如大数定律等,都以这个为前提。
- 什么叫随机实验?随机实验需要满足3个条件:①相同条件下,可重复进行; ②一次试验,结果不确定,但所有结果已知(比如抛硬币,要么正要么反,所有结果已知); ③试验之前,结果未知。
2. 先验、似然、后验概率
大家可能都清楚,就是先验是因,后验是果,似然就是在因的条件下,发生果的概率。那么到底什么是因,什么是果呢?
还是举个例子讲一下先验、似然、观测、后验。
比如在相机和IMU融合SLAM中(由于博主是做这个的,就以这个为例子了,可能对不是做这个的不友好,所以博主上面给了参考博客链接),机器人的状态(包括位姿、速度、偏置等)的概率分布是先验,相机拍摄的图像信息(特征点、特征匹配等)和IMU提供的运动信息(加速度、角速度积分成速度、位移、旋转),都是观测,你得先有机器人的状态,才会有这些观测。
P(状态) - 先验 , 记为 P(A)
P(状态/图像和运动信息) - 后验,在知道观测的条件下,状态的概率(知果求因) , 记为P(A/B)
P(图像和运动信息/状态) - 似然概率,在知道状态的条件下,发生这些观测的概率(知因求果),记为P(B/A)
贝叶斯公式: P(B/A)∗P(A)=P(A/B)∗P(B), 变形为:P(B/A)=P(A)(P(A/B)∗P(B)
若有多个观测,那么P(Bj/A)=∑i=1i=nP(A/Bi)∗P(Bi)(P(A/Bj)∗P(Bj)
有一个问题,不知道大家是否想过,为什么状态会是一个概率,难道它是随机事件吗?
个人认为,是因为我们通常会认为我们估计出来的状态,由于传感器观测会包含随机噪声,所以状态也可认为服从一个高斯分布,它的均值就是最优的状态结果。
二、贝叶斯滤波
1. 问题建模
(1) 符号定义:
- X 表示先验,也就是状态变量,x 表示状态变量的一个取值。Y 表示观测,y 表示观测的一个取值。
- Q 表示建模误差,R表示传感器误差,都认为服从高斯分布。
- fk(x)表示第k时刻的状态概率密度函数(pdf),它是P(Xk<x)对x的导数。
fk(x)=dxdP(Xk<x)
- fQk(x)表示预测噪声Qk的概率密度函数,fRk(x)表示观测噪声(传感器精度)Rk的概率密度函数
- 下面推导中,注意fk(x)、f(x)和fYk∣Xk(yk∣xk)的区别,f(x)里面的f是状态转移函数,就是下面状态方程里面那个f ; 而fk(x)表示k时刻状态量Xk概率密度函数,同理,fk−1(x)就表示Xk−1概率密度函数 ; 如果是fYk∣Xk(yk∣xk),就表示P(Yk∣Xk)的概率密度函数。 它们没有关系。
(2) 问题描述
- 问题一:我们求的是状态,为什么扯这么多概率?
因为状态服从一个高斯分布,没有贝叶斯滤波前,随机性很高噪声很大,通过这些预测和观测我们把噪声很大的状态的先验概率分布,转化为经过预测和观测纠正的后验概率,降低它的随机性和噪声,噪声的降低体现在概率分布中的方差减小了(后验概率方差<先验概率方差)。后验概率分布的均值就是我们要求的最优状态估计结果。
Xk表示状态xk的概率分布,fk(x)表示状态xk的概率密度函数,我们要求k时刻状态xk的最优估计,那就是求Xk的均值即可。那
μk=∫−∞+∞x∗fk(x)dx
如果fk(x)是高斯分布,那可以从fk(x)的公式中直接获取。
所以我们的求解步骤: 求解fk+(x) -> 计算均值得到最优xk
- 问题二:怎么求解fk+(x)呢?
现在我们有的数据是这样的:k-1时刻的状态xk−1已知 ; 对xk有传感器的观测数据
两部分信息我们都不想浪费掉,因此我们分两步计算fk+(x)
1] 根据上一时刻状态预测,称为预测步,对应着状态方程,预测得到fk−(x)。
2] 根据当前时刻的观测来估计,称为更新步,对应着观测方程,更新得到fk+(x)。
(3)状态方程和观测方程
{Xk=f(Xk−1)+QkYk=h(Xk)+Rk
假设X0、Q1、Q2...Qk、R1、R2...Rk相互独立。
状态量可能是一个向量,因为状态量通常不只有一个,在机器人里面,就包含位姿、速度、偏置等,因此大多数情况下,状态量是一个向量,那么Qk、Rk都是矩阵。另外,如果有多个传感器,那么每个时刻都会有多个观测,这里只写了一个传感器时候的例子。
现在有一组观测值,y1、y2...yk分别表示时刻x1、x2...xk对应的观测(或者叫传感器测量数据,R其实是表示的是传感器的精度)。
2、预测步推导(求先验)
我们首先要根据上一时刻k-1的状态来预测当前状态,得到fk−(x),其实就是求先验
fk(x)=dxdP(Xk<x)
要求解fk(x),可以先求解P(Xk<x),求解方法类似与泰勒一阶展开,展开到上一时刻状态量。
P(Xk<x)P(Xk=u)所以,P(Xk<x)因此:fk−(x)=u=−∞∑xP(Xk=u),u连续取值=v=∞∑+∞P(Xk=u∣Xk−1=v)P(Xk−1=v)=v=∞∑+∞P[Xk−f(Xk−1)=u−f(v)∣Xk−1=v]∗P(Xk−1=v)=v=∞∑+∞P[QK=u−f(v)∣Xk−1=v]∗P(Xk−1=v)=v=∞∑+∞P[QK=u−f(v)]∗P(Xk−1=v),因为Xk−1和Qk独立=ε−>0limv=∞∑+∞fQk[u−f(v)]∗ε∗fk−1(v)∗ε,见下面注解1=ε−>0lim∫−∞+∞fQk[u−f(v)]∗fk−1(v)dv∗ε,见下面注解2=u=−∞∑xε−>0lim∫−∞+∞fQk−1[u−f(v)]fk−1(v)dv∗ε=∫−∞x∫−∞+∞fQk[u−f(v)]∗fk−1(v)dvdu=dxdP(Xk<x)=∫−∞+∞fQk[x−f(v)]∗fk−1(v)dv
-
注解1,这里顺便把连续贝叶斯公式推导一下,看完推导也就明白上面公式了:
P(X<x∣Y=y)又因为:P(X<x∣Y=y)因此:∫−∞xfX∣Y(x∣y)dx所以,fX∣Y(x∣y)=u=−∞∑x(P(X=u∣Y=y)=u=−∞∑xp(Y=y)P(Y=y∣X=u)P(x=u)=ε−>0limu=−∞∑xP(y<Y<y+ε)P(y<Y<y+ε∣X=u)P(u<X<u+ε)=ε−>0limu=−∞∑xfY(η3)∗εfY∣X(η1∣u)∗ε∗fX(η2)∗ε,其中η1∈(y,y+ε),η2∈(u,u+ε),η3∈(y,y+ε)=ε−>0limu=−∞∑xfY(y)fY∣X(y∣u)∗fX(u)∗ε,因为η1、η3∈(y,y+ε),而ε→0,所以取η1、η3=y,η2取u=∫−∞xfY(y)fY∣X(y∣u)∗fX(u)du,这里利用了定积分定义,参考注解2的知乎链接=∫−∞xfY(y)fY∣X(y∣x)∗fX(x)dx=∫−∞xfX∣Y(x∣y)dx=∫−∞xfY(y)fY∣X(y∣x)∗fX(x)dx对连续随机变量,也有:=fY(y)fY∣X(y∣x)∗fX(x)(贝叶斯公式)
-
注解2:参考链接:利用定积分定义求极限的原理与套路,你会了吗?
3、更新步推导(求后验)
假设有一个传感器,在k时刻有观测Yk=yk,我们需要利用观测去更新fk−(x)成fk+(x),其实就是后验概率fk(x∣yk)
fYk∣Xk(yk∣x)所以,有:fk+(x)其中,η为:η=ε→0limεP(yk<Yk<yk+ε∣Xk=x)=ε→0limεP(yk−h(x)<Yk−h(Xk)<yk−h(x)+ε∣Xk=x)=ε→0limεP(yk−h(x)<Rk<yk−h(x)+ε∣Xk=x)=ε→0limεP(yk−h(x)<Rk<yk−h(x)+ε),因为Rk和Xk独立=fRk[yk−h(x)]=fk(x∣yk)=fYk(yk)fYk∣Xk(yk∣x)∗fk−(x)=fYk(yk)fRk[yk−h(x)]∗fk−(x)=η∗fRk[yk−h(x)]∗fk−(x)={∫−∞+∞fRk([yk−h(x)]∗fk−(x)dx}−1
经过更新步之后,状态量的方差大大降低,也就是经过贝叶斯滤波后减少了噪声,得到更优更可靠的状态估计值.
三、贝叶斯滤波递推过程
f0(x)其中,ηk预测f1−(x)=∫−∞+∞fQ1[x−f(v)]f0(v)dv观测更新f1+(x)=η1∗fR1[y1−h(x)]∗f1−(x)预测f2−(x)=∫−∞+∞fQ2[x−f(v)]f1+(x)dv观测更新f2+(x)=η2∗fR2[y2−h(x)]∗f2−(x)……预测fk−(x)=∫−∞+∞fQk[x−f(v)]fk−1+(x)dv观测更新fk+(x)=ηk∗fRk[yk−h(x)]∗fk−(x)={∫−∞+∞fRk[yk−h(x)]∗fk−(x)dx}−1,其实就是后面那一坨的积分的倒数。
公式中,除了初值即可认为是先验,又可以认为是后验,其他都需要经过预测和更新两个步骤,实验证明,初值即使很差,也不会影响后面时刻状态量的估计值。
这里递推的都是概率密度函数,要求各时刻的状态估计,只需要求均值即可。
x^k=∫−∞+∞x∗fk+(x)dx
四、贝叶斯滤波算法流程
(1) 设初值:初始状态x0和它的概率密度函数f0(x)
(2) 预测步:fk−(x)=∫−∞+∞fQk[x−f(v)]∗fk−1+(v)dv
(3) 更新步:fk+(x)=ηk∗fRk[yk−h(x)]∗fk−(x),其中,ηk={∫−∞+∞fRk[yk−h(x)]∗fk−(x)dx}−1
这里要提一下有些博客不准确的说法:P(xk∣yk)=η∗P(yk∣xk)∗∫−∞+∞P(xk∣xk−1)∗P(xk−1)dxk−1
(4) 求状态量:x^k=∫−∞+∞x∗fk+(x)dx
五、贝叶斯滤波算法优缺点
- 可以有效滤除噪声,得到比较精准的状态估计;
- 缺点也很明显,从fk−1+(x)→fk−(x),计算η,计算期望x^k,都需要做无穷积分,大多数情况下没有解析解。
到这里,贝叶斯滤波就讲完了,后面会写一写关于卡尔曼滤波的博客,搞清楚贝叶斯滤波之后,卡尔曼滤波就非常好理解了,因为它就是在贝叶斯滤波基础上,将状态方程里面的函数f和观测方程里的h假设是线性的,没有了非线性,计算也就发便多了。
标签:Xk,yk,infty,fY,滤波,贝叶斯,dx,搞懂,fk
来源: https://blog.csdn.net/wq1psa78/article/details/105849353