非参数滤波是高斯滤波的一种替代选择,非参数滤波不依赖于后验,它是通过有限的数据来近似后验。近似的质量主要取决于表示后验数据的量,当数据量足够大的时候,往往会一致收敛于真实的后验值。本章主要讨论了两种非参数的滤波,直方图滤波和粒子滤波。

直方图滤波将状态空间分解为有限多个区域,并且用一个概率值表示每一个区域的累积后验。当应用于有限空间时,被称为离散贝叶斯滤波,应用于连续空间时被称为直方图滤波。机器人技术中有许多的分解技术,可以分为依赖和不依赖环境结构这两种分解技术,当分解依赖于环境结构时,称为拓扑。也可分为静态和动态的分解技术。静态分解需要提前进行,而不需要考虑置信度。动态分解依赖于置信度的特性,经常按照后验概率成比例的提高空间分辨率。动态分解能够给出较好的分辨结果,但是实现动态分解比静态分解更难。

粒子滤波的主要思想是用一系列从上一状态后验得到的随机状态中采样得到新的后验。这种算法具有比较好的鲁棒性,而且可以表示比高斯分布更为广泛的分布空间。有一些策略可以减少粒子滤波的误差,其中最流行的方法是减少由算法的随机性引起的估计方差技术和根据后验复杂性调节粒子数量的技术。

离散贝叶斯滤波算法

离散贝叶斯滤波算法其实就是将积分改为求和,思路也非常的简单。直接给出算法过程吧。

整个算法如下:

对于连续空间而言就是将整个空间分成若干个直方图区域,将同一个直方图区域的概率看成均匀分布的。

通过这个思路,我们就可以写出直方图滤波的算法。

静态二值贝叶斯滤波

机器人技术中某些问题可以表达为不随时间变化的二值状态的最优估计问题。如果一个机器人从传感器测量的序列中估计一个环境是一个二值问题,比如估计环境中的门是开还是闭这样的问题。

当状态静止的时候,置信度就是测量函数,状态只可能存在两个值和,两者之和为1。这种情况可以用离散贝叶斯来计算,但是也可以用概率比对数来实现。状态x的概率定义为此事件发生的概率除以此事件不发生的概率:

概率对数就是取对数:

概率对数假设值为整个实数,从负无穷到正无穷,更新以概率对数表示的置信度的贝叶斯滤波计算很简单!这里通过使用反向测量模型p(x|zt) 代替了熟悉的前向模型 p(zt|x) 。这里通俗的解释一下,什么是反向观测模型。比如从相机图像中观测门的开关情况,这里状态只有两个,但是需要测量的整个空间却非常大,我们通过设计一个函数根据相机图像来计算门关着的概率比去计算门关着的概率分布会更简单。也就是说我们实现一个反向测量比实现前向测量要容易很多。

整个算法如下:

非常简单,就一个式子。接下来是重点中的重点了,粒子滤波

粒子滤波

粒子滤波器是贝叶斯滤波的一种非参数实现,它也是通过有线个参数来近似后验,但是这些参数的生成方式与直方图不一样。粒子滤波器用一些列来自该分布的样本来表示这个分布。再粒子滤波器中,后验分布的样本叫做粒子滤波:

也就是说状态空间的一个子区域被样本填充得越多,那么真实状态位于该区域的可能性就越大。

先给出一个粒子滤波的基本算法思路:

重采样是一种基于达尔文适者生存思想的概率实现:它将粒子集重新聚集在状态空间中具有较高后验概率的区域,也就是说它把计算资源较为集中的分配给了概率高的区域。

粒子滤波的误差主要来源于随机采样。因为从有限的样本去估计概率密度始终是会存在误差的,也就是说样本的统计特性和初始密度的统计特性是存在偏差的,因此会造成误差,这样产生的误差也叫抽样方差。最开始的时候,粒子集由先验信息产生且粒子集分布在整个状态空间,但是重采样有时候并不能完全复现。随着时间的推移,越来越多的粒子由于重采样被简单随机的剔除,并没有加入新的状态(因为粒子是不断的聚集在概率较大的状态空间里面,所以这里并没有引入新的状态),这种情况下最极端的可能就是仅剩一种状态,概率就变为了1。也就是说尽管粒子本身的不一致性减小了,但是作为真实置信度估计器的粒子群不一致性却增加了(因为初始粒子有很多状态,到了后面状态在不断的被剔除)。因此控制粒子滤波器的这种变化或者误差是非常有必要的。

目前有两种方法减少粒子群与真实置信度估计器之间的差异。第一种是减少采样频率,重采样太频繁可能导致丢失多样性,重采样频率太低则会浪费过多的时间在低概率区域。如何决定重采样通常是测量权值的变化,如果权值变化为0,那就不需要进行重采样(因为没有发生变化),如果权值变化较大则需要进行重采样。

减少误差的第二种策略叫做低方差重采样。该算法并不是选择M个随机数以及他们对应的粒子,而是计算一个单独的随机数并且根据这个随机数选择样本,同时也必须满足概率和权重成正比.。

下面,以低方差重采样给出粒子滤波的matlab算法程序:

clc; 
M=1000;%粒子数 
P0=5;  %初始状态协方差 
Q=10;   %过程噪声方差 
R=4;    %量测方差阵
tf=200; %终止时间  
pdf_v=inline('1/(2*pi*1)^(1/2)*exp(-(x.^2)/(2*1))'); 
f=inline('x./2+25*x./(1+x.^2)+8*cos(1.2*t)','x','t');%状态转移方程
h=inline('(x.^2)/20');       %量测方程
x(1)=sqrtm(P0)*randn(1);       %初始状态值,sqrtm矩阵开根函数,randn(n)返回月一个n*n的随机矩阵
y(1)=feval(h,x(1))+sqrtm(R)*randn(1);  
for t=2:tf    %系统仿真    
    x(t)=feval(f,x(t-1),t-1)+sqrtm(Q)*randn(1);  
    y(t)=feval(h,x(t))+sqrtm(R)*randn(1); 
end
xTrue=x;  
xhat=PF(f,h,pdf_v,Q,P0,M,y);%状态值、量测值、高斯分布、过程噪声方差、初始方差阵、粒子数、包含噪声的量测值 
plot(1:tf,xhat,'b--',1:tf,xTrue,'r'); 
xlabel('时间');  
legend('状态估计值','状态真实值');
title('粒子滤波仿真实验'); 
grid on;

function xhat=PF(f,h,pdf_v,Q,P0,M,y) %状态值、量测值、高斯分布、过程噪声方差、
                                     %初始方差阵、粒子数、包含噪声的量测值
n=size(P0,2); %返回P0的列数
x=sqrtm(P0)*randn(n,M);%初始化粒子
tf=size(y,2); 
for t=1:tf    
    e=repmat(y(t),1,M)-h(x);        %计算权重   
    w=feval(pdf_v,e);               %转为概率    
    w=w/sum(w);   %归一化权值 
   
    xhat(t)=sum(repmat(w,n,1).*x,2);%行求和   
    ind=resampling(x,w);              %重采样  
  
    x=ind;%新粒子
    x=feval(f,x,t)+sqrtm(Q)*randn(n,M);%时间更新 
end
end
function [i]=resampling(x,w)

wc=w(1);
M=length(w); 
r=rand*(1/M);
%u=rand(1,M)/M;
i=zeros(1,M);
k=1;
for j=1:M 
    u=r+(j-1)/M;
    while(wc<u)      
        k=k+1;
        wc=wc+w(k);
    end
    i(j)=x(k); 
end
end

源代码来自百度文库粒子滤波程序一看就懂_百度文库,略有改动

最终滤波效果如下: