2 粒子滤波理论
粒子滤波通过非参数化的蒙特卡洛(Monte Carlo)模拟方法来实现递推贝叶斯滤波，适用于任何能用状态空间模型描述的非线性系统，精度可以逼近最优估计。粒子滤波器具有简单、易于实现等特点，它为分析非线性动态系统提供了一种有效的解决方法，从而引起目标跟踪、信号处理以及自动控制等领域的广泛关注。本章首先概述用于求解目标状态后验概率的贝叶斯滤波理论，随后介绍具有普遍适用性的粒子滤波器，最后针对当前粒子滤波器存在的粒子多样性丧失问题，提出了一种量子进化粒子滤波算法。
2.1 贝叶斯滤波
动态系统的目标跟踪问题可以通过图2.1所示的状态空间模型来描述。本节在贝叶斯滤波框架下讨论目标跟踪问题。
图2.1 状态空间模型 Fig. 2.1 State space model
在目标跟踪问题中，动态系统的状态空间模型可描述为
xk?f(xk?1)?uk?1yk?h(xk)?vk (2.1)
其中f(?),h(?)分别为状态转移方程与观测方程，xk为系统状态，yk为观测值，uk为过程
vk为观测噪声。噪声，为了描述方便，用Xk?x0:k?{x0,x1,?,xk}与Yk?y1:k?{y1,?,yk}分别表示0到k时刻所有的状态与观测值。在处理目标跟踪问题时，通常假设目标的状态转移过程服从一阶马尔可夫模型，即当前时刻的状态xk只与上一时刻的状态xk-1有关。另外一个假设为观测值相互独立，即观测值yk只与k时刻的状态xk有关。
贝叶斯滤波为非线性系统的状态估计问题提供了一种基于概率分布形式的解决方案。贝叶斯滤波将状态估计视为一个概率推理过程，即将目标状态的估计问题转换为利用贝叶斯公式求解后验概率密度p(Xk|Yk)或滤波概率密度p(xk|Yk)，进而获得目标状态的最优估计。贝叶斯滤波包含预测和更新两个阶段，预测过程利用系统模型预测状态的先验概率密度，更新过程则利用最新的测量值对先验概率密度进行修正，得到后验概率密度。
假设已知k?1时刻的概率密度函数为p(xk?1|Yk?1)，贝叶斯滤波的具体过程如下： (1) 预测过程，由p(xk?1|Yk?1)得到p(xk|Yk?1)：
p(xk,xk?1|Yk?1)?p(xk|xk?1,Yk?1)p(xk?1|Yk?1) (2.2)
当给定xk?1时，状态xk与Yk?1相互独立，因此
p(xk,xk?1|Yk?1)?p(xk|xk?1)p(xk?1|Yk?1) (2.3)
上式两端对xk?1积分，可得Chapman-Komolgorov方程
p(xk|Yk?1)??p(xk|xk?1)p(xk?1|Yk?1)dxk?1
(2.4)
(2) 更新过程，由p(xk|Yk?1)得到p(xk|Yk)：
获取k时刻的测量yk后，利用贝叶斯公式对先验概率密度进行更新，得到后验概率
p(yk|xk,Yk?1)p(xk|Yk?1)p(yk|Yk?1) p(xk|Yk)? (2.5)
假设yk只由xk决定，即
p(yk|xk,Yk?1)?p(yk|xk) (2.6)
因此
p(xk|Yk)?p(yk|xk)p(xk|Yk?1)p(yk|Yk?1) (2.7)
其中，p(yk|Yk?1)为归一化常数
p(yk|Yk?1)??p(yk|xk)p(xk|Yk?1)dxk (2.8)
贝叶斯滤波以递推的形式给出后验(或滤波)概率密度函数的最优解。目标状态的最优估计值可由后验(或滤波)概率密度函数进行计算。通常根据极大后验(MAP)准则或最小均方误差(MMSE)准则，将具有极大后验概率密度的状态或条件均值作为系统状态的估计值，即
?kMAPx =argminp(xk|Yk) (2.9)
xkSE x ) d (2.10) ?kM Mx =fE[x(Y)?|?]fx(pkx)(kY| kkkk贝叶斯滤波需要进行积分运算，除了一些特殊的系统模型（如线性高斯系统，有限状态
的离散系统）之外，对于一般的非线性、非高斯系统，贝叶斯滤波很难得到后验概率的封闭解析式。因此，现有的非线性滤波器多采用近似的计算方法解决积分问题，以此来获取估计的次优解。在系统的非线性模型可由在当前状态展开的线性模型有限近似的前提下，基于一阶或二阶Taylor级数展开的扩展Kalman滤波得到广泛应用[119]。在一般情况下，逼近概率密度函数比逼近非线性函数容易实现。据此，Julier与Uhlmann提出一种Unscented Kalman滤波器，通过选定的sigma点来精确估计随机变量经非线性变换后的均值和方差，从而更好的近似状态的概率密度函数，其理论估计精度优于扩展Kalman滤波一中方案便是基于蒙特卡洛模拟的粒子滤波器。
[120]
。获取次优解的另外
2.2 粒子滤波
早在20世纪50年代，Hammersley便采用基于序贯重要性采样(Sequential importance sampling,SIS)的蒙特卡洛方法解决统计学问题[121]。20世纪60年代后期，Handschin与Mayne使用序贯蒙特卡洛方法解决自动控制领域的相关问题[122]。20世纪70年代，Handschin、Akashi以及Zaritskii等学者的一系列研究工作使得序贯蒙特卡洛方法得到进一步发展 [126]
。限于当时的计算能力以及算法本身存在的权值退化问题，序贯重要性采样算法没有受到足够重视，在随后较长一段时间内进展较为缓慢。直到20世纪80年代末，计算机处理能力的巨大进展使得序贯蒙特卡洛方法重新受到关注。Tanizaki、Geweke等采用基于重要性采样的蒙特卡洛方法成功解决了一系列高维积分问题[127-130]。Smith与Gelfand提出的采样-重采样思想为Bayesian推理提供了一种易于实现的计算策略[131]。随后，Smith与Gordon等人合作，于20世纪90年代初将重采样(Resampling)步骤引入到粒子滤波中，在一定程度上解决了序贯重要性采样的权值退化问题，并由此产生了第一个可实现的SIR(Sampling importance resampling)粒子滤波算法(Bootstrap滤波)[132]，从而掀起粒子滤波的研究热潮。美国海军集成水下监控系统中的Nodestar便是粒子滤波应用的一个实例。进入21世纪，粒子滤波器成为一个非常活跃的研究领域，Doucet、Liu、Arulampalam等对粒子滤波的研究作了精彩的总结，IEEE出版的论文集“Sequential Monte Carlo Methods in Practice”对粒子滤波器进行了详细介绍[136]。
[133-135]
[123][124, 125]
2.2.1 贝叶斯重要性采样
蒙特卡洛模拟是一种利用随机数求解物理和数学问题的计算方法，又称为计算机随机模拟方法。该方法源于第一次世界大战期间美国研制原子弹的曼哈顿计划，著名数学家冯?诺伊曼作为该计划的主持人之一，用驰名世界的赌城，摩纳哥的蒙特卡洛来命名这种方法。蒙特卡洛模拟方法利用所求状态空间中大量的样本点来近似逼近待估计变量的后验概率分布，如图2.2所示，从而将积分问题转换为有限样本点的求和问题。粒子滤波算法的核心思想便是利用一系列随机样本的加权和表示后验概率密度，通过求和来近似积分操作。假设可以从
(i)后验概率密度p(xk|Yk)中抽取N个独立同分布的随机样本xk，i?1,?,N，则有
p(xk|Yk)?1NN??(xi?1k(i)?xk) (2.11)
这里xk为连续变量，?(x-xk)为单位冲激函数(狄拉克函数)，即?(x-xk)?0,x?xk，且??(x)dx?1。当xk为离散变量时，后验概率分布P(xk|Yk)可近似逼近为
1NNP(xk|Yk)???(xi?1k(i)?xk) (2.12)
其中，?(xk?xk(i))?1,xk?xk(i); ?(xk?xk(i))?0,xk?xk(i)。
P(x)1
图2.2 经验概率分布函数
Fig. 2.2 Empirical probability distribution function
(i)设xk为从后验概率密度函数p(xk|Yk)中获取的采样粒子，则任意函数f(xk)的期望
估计可以用求和方式逼近，即
E[f(xk)|Yk]??f(xk)p(xk|Yk)dxk?1NN?i?1f(xk) (2.13)
(i)蒙特卡洛方法一般可以归纳为以下三个步骤：
(1)构造概率模型。对于本身具有随机性质的问题，主要工作是正确地描述和模拟这个概率过程。对于确定性问题，比如计算定积分、求解线性方程组、偏微分方程等问题，采用蒙特卡洛方法求解需要事先构造一个人为的概率过程，将它的某些参量视为问题的解。
(2)从指定概率分布中采样。产生服从己知概率分布的随机变量是实现蒙特卡洛方法模拟试验的关键步骤。 (3)建立各种估计量的估计。一般说来，构造出概率模型并能从中抽样后，便可进行现模拟试验。随后，就要确定一个随机变量，将其作为待求解问题的解进行估计。 在实际计算中，通常无法直接从后验概率分布中采样，如何得到服从后验概率分布的随机样本是蒙特卡洛方法中基本的问题之一。重要性采样法引入一个已知的、容易采样的重要性概率密度函数q(xk|Yk)，从中生成采样粒子，利用这些随机样本的加权和来逼近后验滤
波概率密度p(xk|Yk)，如图2.3所示。令{xk(i),wk(i),i?1,?.N}表示一支撑点集，其中xk(i)为是k时刻第i个粒子的状态，其相应的权值为wk(i)，则后验滤波概率密度可以表示为
N p(xk|Yk)?其中，
?wi?1(i)k?x(k?xki() ) (2.14)
w??(i)kp(xk|Yk)q(x(i)k(i)|Yk) (2.15)
图2.3 重要性采样 Fig. 2.3 Importance sampling
当采样粒子的数目很大时，式(2.14)便可近似逼近真实的后验概率密度函数。任意函数
f(xk)的期望估计为
E[f(xk)|Yk]=1NN?i?1f(x(i)k)p(xk|Yk)q(xk|Yk)(i)(i)?1NN?i?1f(xk)wk (2.16)
(i)(i)2.2.2 序贯重要性采样算法
在基于重要性采样的蒙特卡洛模拟方法中，估计后验滤波概率需要利用所有的观测数据，每次新的观测数据来到都需要重新计算整个状态序列的重要性权值。序贯重要性采样作为粒子滤波的基础，它将统计学中的序贯分析方法应用到的蒙特卡洛方法中，从而实现后验滤波概率密度的递推估计。假设重要性概率密度函数q(x0:k|y1:k)可以分解为
q(x0:k|y1:k)?q(x0:k?1|y1:k?1)q(xk|x0:k?1,y1:k) (2.17)
设系统状态是一个马尔可夫过程，且给定系统状态下各次观测独立，则有
kp(x0:k)?p(x0)?p(xi|xi?1) (2.18)
i?1