《第二章-贝叶斯状态估计与粒子滤波(共30页).doc》由会员分享，可在线阅读，更多相关《第二章-贝叶斯状态估计与粒子滤波(共30页).doc（31页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上第二章 贝叶斯状态估计与粒子滤波视觉跟踪可视为状态估计问题16,54，即根据视觉目标在先前帧的状态信息估计其在当前帧的状态，从而实现视觉跟踪。状态估计一直都是自动控制、通讯、航空与航天等领域的经典研究主题之一69,70。贝叶斯状态估计是处理不确定性条件下状态估计问题的有力理论工具21,22,71。为了有效处理非高斯、非线性状态估计问题，二十世纪九十年代人们提出了粒子滤波19-22,71，粒子滤波是基于Monte Carlo随机模拟的贝叶斯滤波方法。本章将简单介绍贝叶斯状态估计和粒子滤波相关理论问题。首先，通过介绍贝叶斯状态估计相关理论，引出贝叶斯状态滤波问题及实现贝叶
2、斯状态滤波的两大理论工具：卡尔曼系滤波器和粒子滤波。然后，简单介绍了卡尔曼系滤波器的相关理论和算法。最后，详细介绍了粒子滤波理论框架、收敛性问题及经典采样策略。2.1 贝叶斯状态估计估计理论是概率论和数理统计的一个分支，所研究的对象是随机现象。它是根据受干扰的观测数据来估计关于随机变量、随机过程或系统的某些特性的一种数学方法70。所谓估计，就是从带随机噪声干扰的观测信号中提取有用信息，可定义如下：定义 2.1 如果假设被估计量为维向量，而其观测量为维向量，且观测量与被估计量之间具有如下关系 (2.1)其中，是已知的维向量函数，由观测方法决定；是观测误差向量，通常为一个随机过程。那么，所谓估计问
3、题，就是在时间区间内对进行观测，从而在得到观测数据的情况下，要求构造一个观测数据的函数去估计的问题，并称是的一个估计量，或称的估计为69,70。一般地，估计问题可以分为两类：状态估计和参数估计。状态和参数的基本差别在于，前者是随时间变化的随机过程，后者是不随时间变化或随时间缓慢变化的随机变量。因此，可以说状态估计是动态估计，而参数估计是静态估计。在此，主要讨论系统状态估计问题。下面首先介绍最优估计和估计准则，然后在贝叶斯意义下描述状态估计问题。2.1.1 最优估计与估计准则一般地，在实际工程应用中总希望估计出来的参数或状态愈接近真实值愈好，即如何最优地利用系统观测数据得到一个最优估计量，这就是
4、最优估计问题。所谓最优状态估计，是指在某一确定的估计准则条件下，按照某种统计意义，使系统状态估计达到最优70。因此，最优状态估计是针对某一估计准则而言的。估计准则是衡量估计的好坏的，选择合理的估计准则是极其重要的。可以说，估计准则在很大程度上决定了估计的性能、求解估计问题所使用的估计方法及估计量的性质等。估计准则是多种多样的，但在贝叶斯意义下，统计估计准则可分为：贝叶斯估计准则与非贝叶斯估计准则69,70。非贝叶斯估计准则常见的有最小二乘准则和极大似然准则；而贝叶斯估计准则常见的有极大后验准则和最小方差准则。估计准则和估计方法是紧密相关的，选择不同的估计准则就对应不同的估计方法。下面将简要介绍
5、基于这几种估计准则的统计估计方法：1. 最小二乘估计最小二乘估计是法国数学家高斯（Gauss）于1809年提出的，是一种使用最广泛的估计方法之一。最小二乘估计可定义如下：定义 2.2 设被估计量是非时变的维随机向量，如果对其进行次线性观测，则有 (2.2)其中是维观测向量，是观测矩阵，是维的零均值观测误差向量。如果将次线性观测简写成， (2.3)其中，则是一个维向量，是矩阵，是维的零均值观测误差向量。当时，可根据来估计。如果要求选择的一个估计，使得性能指标 (2.4)或更一般形式的二次型性能指标 (2.5)达到极小，那么称这个估计为的最小二乘估计或加权最小二乘估计，并记为或。其中为对称正定加权
6、矩阵69,70。对于最小二乘估计，作如下几点说明：（1）是所有观测数据的全体，因此最小二乘估计要求把所有观测数据都储存起来作统一处理，很难实现实时处理。（2）最小二乘估计或加权最小二乘估计都是无偏估计。（3）设观测误差的方差阵为，则可以证明，当选择加权矩阵时，能使加权最小二乘估计的方差阵达到最小70。2. 极大似然估计极大似然估计是以观测值出现的概率最大作为准则的，是应用非常广泛的参数估计方法。1906年费希尔（Fisher）首先使用这种估计方法，它是以似然函数概念为基础的。极大似然估计可定义如下69,70,72,73：定义 2.3 设为维被估计量，为的次观测数据集，它是从同一个分布独立采样得
7、到的（即独立同分布的）。记 (2.6)称为样本的似然函数。如果样本集的一个函数满足：，则成为的极大似然估计70。对于极大似然估计，可以证明，当观测次数趋于无限大时，极大似然估计量是一种无偏估计量。因此，极大似然估计是渐近无偏的。此外，极大似然估计量可以是随机量，也可以是非随机参数，适用范围较广。3. 贝叶斯估计贝叶斯估计是以贝叶斯统计为基础的，是当前最优估计的研究热点之一。贝叶斯统计的主要优势在于能处理数据分析中的不确定性，包括被估计的参数以及模型中的任何不确定性。在贝叶斯统计中，被估计量都被当作随机变量，它服从一定的分布，并认为已观测到的数据可以揭示这个分布的信息。根据先验知识确定的被估计量
8、的分布，称为先验分布；该分布将根据纳入的观测数据中的信息进行改进而得到后验分布。从先验分布进化为后验分布是通过贝叶斯定理来实现的。根据贝叶斯统计原理，贝叶斯估计可定义如下：定义 2.4 设为被估计量，其先验分布为，为的个观测值，其条件概率密度函数为，则可利用贝叶斯公式得到后验概率密度函数 (2.7)如果令为损失函数，度量的一个估计带来的损失，则后验分布下的损失函数的期望 (2.8)称为贝叶斯后验风险。如果估计能使后验风险达到最小，则称为的Bayes估计70。由式(2.8)可知，损失函数的选择非常重要，选择不同形式的损失函数，就可得到不同的贝叶斯估计结果。下面将讨论两种情况：最大后验估计和最小方
9、差估计。（1）最大后验估计：如果令损失函数为风险，即 (2.9)于是，将式(2.9)代入式(2.8)则有 (2.10)显然，由式(2.10)可知，使贝叶斯后验风险达到最小，就相当于要求后验概率密度达到最大70。也就是说，“达到最小”与“达到最大”是等价的。因此，可以把“使后验概率密度函数达到最大”作为估计准则来得到贝叶斯估计，并把这种估计称为最大后验估计，记为。对于最大后验估计，将式(2.7)代入式(2.10)，即有 (2.11)由式(2.11)可知，在对没有任何先验统计知识的情况下，最大后验估计就退化为极大似然估计，因此可以说，极大似然估计是一种特殊的最大后验估计，或者说是一种特殊的贝叶斯估
10、计。但是，在一般情况下，由于考虑了的先验知识（即先验分布），因此最大后验估计将优于极大似然估计。（2）最小方差估计：如果令估计误差的二次型函数为损失函数，即 (2.12)于是，将式(2.12)代入式(2.8)则可得到贝叶斯性能指标为 (2.13)一般地，把“使式(2.13)所示的贝叶斯性能指标达到最小”作为估计准则所得到的贝叶斯估计称为最小方差估计，记为。对于最小方差估计，作如下几点说明：、最小方差估计为无偏估计。、可以证明，任何其他估计的均方误差阵或任何其他无偏估计的方差阵都将大于最小方差估计的误差方差阵，即最优估计具有最小的估计误差方差阵70。2.1.2 贝叶斯意义下的状态估计卡尔曼（Ka
11、lman）开创性地将状态变量和状态空间概念引入到最优估计，提出了状态估计理论74。从状态空间观点，状态是比信号更为广泛、更灵活的概念，非常适合处理多变量系统，信号可以视为状态或状态的分量。一般地，动态系统的状态可定义如下：定义 2.5 把能完全确定动态系统运动特性的最小一组变量，称为该动态系统的状态；状态变量所构成的向量称为状态向量；状态向量所张成的空间称为状态空间。对于实际系统而言，其状态很难直接获得或不允许直接测量，得到的只是与状态有关的一些观测数据，而且观测数据往往会受到随机噪声的干扰，是有观测误差的。为了得到系统的状态，就只有根据这些观测数据构造或估计系统的状态，当然，系统状态的估计值
12、应尽量接近实际状态，这就是所谓的系统状态估计问题69,70。在统计理论里，贝叶斯统计是处理不确定性问题的有力工具75。一直以来，贝叶斯状态估计是系统状态估计（特别是非线性、非高斯状态估计）的一大研究热点21,48,71,76。下面将简要介绍离散系统的状态估计和贝叶斯意义下的状态估计：1. 离散系统的状态估计离散系统的状态估计可定义如下：定义 2.6 设随机、离散系统S的状态空间模型为 (2.14)其中，为时刻（）的系统状态向量，为系统随机噪声，为系统状态转移模型；为时刻的系统观测向量，为随机观测噪声，为系统观测模型。如果对系统的状态向量进行次观测，从而得到观测序列，那么所谓离散系统得状态估计问
13、题，就是要求根据整个观测数据，求得在时刻系统状态向量的一个最优估计量的问题，通常把所得到的估计量记为，并且按照和的不同关系，状态估计可分为三类：1） 称为预测（或外推）；2） 称为滤波；3） 称为平滑（或内插）。69, 702. 贝叶斯意义下的状态估计与递推贝叶斯滤波考虑式（2.14）所示状态空间模型建模的动态系统S的状态估计问题。如果将系统状态转移模型和观测模型概率化，则系统S包含两个随机过程：状态过程和观测过程。其中，状态过程可视为具有初始分布的离散马尔科夫链，且其马尔科夫转移核为，称为状态转移概率；观测过程条件依赖于状态过程，其条件分布为，称为似然函数（又称为似然比）。于是，系统S的状态
14、过程和观测过程的统计特性可完全由初始分布、状态转移概率和似然函数决定： (2.15)对于式（2.14），如果与为独立同分布的噪声，且相互独立，统计特性已知（即，已知），则状态转移概率和似然函数可由噪声的分布给出： (2.16)特别地，当噪声信号具有高斯分布，且均值和协方差阵为：，则式（2.16）可改写为： (2.17)其中，为高斯分布函数。于是，根据贝叶斯估计准则（定义2.4）和系统状态估计（定义2.6），贝叶斯意义下的状态估计问题即为：给定一系列观测数据，求得一列“最优状态估计序列”，以使得式（2.8）定义的贝叶斯后验风险指标达到最小： (2.18)对于许多实际问题，最感兴趣的是状态的当前值
15、。只考虑当前状态取值估计，这样的最优状态估计问题就变成一个最优滤波问题：在给定的观测数据，求得一个最优的当前状态估计值，使得满足贝叶斯后验风险指标： (2.19)由贝叶斯估计（见2.1.1节）可知，对于其两种估计方法（最大后验估计和最小方差估计），若分别知道和的显式解，则贝叶斯意义下的最优估计和最优滤波问题的递推显式解即可求得。一般地，求解贝叶斯估计和滤波的后验分布的递推方程称为递推贝叶斯估计和滤波方程。对于离散时间随机系统S，在给定初始先验状态分布条件下，的递推解可按如下序贯贝叶斯估计公式给出： (2.20)其中，为归一化常数。在实际应用中，感兴趣的是当前状态的滤波估计，如何求得状态后验分布
16、是贝叶斯滤波的核心问题。对于贝叶斯滤波，可由如下两式递推求解：（1）预测方程（Chapman-Kolmogorov方程）： (2.21)（2）更新方程（贝叶斯推理）： (2.22)一般地，式（2.21）和式（2.22）被合称为贝叶斯滤波方程。在理论意义上，式（2.20）与式（2.21）、式（2.22）分别构成了贝叶斯最优估计和贝叶斯最优滤波的递推最优求解方法。但是，在一般情况下，它们的显式解是不能得到的。因此，人们采用各种途径寻找贝叶斯滤波的近似数值求解方法。如果系统具有线性和高斯性，则卡尔曼滤波是求解最优贝叶斯滤波有效方法69,70,74,77,78。为了处理非线性滤波问题，人们提出了扩展卡
17、尔曼滤波69,70,77-80；同时，为了更好地处理非线性滤波问题，Julier等在卡尔曼滤波的理论框架下提出了Unscented卡尔曼滤波52,53。但是，卡尔曼滤波、扩展卡尔曼滤波和Unscented卡尔曼滤波都不能有效解决非线性、非高斯情况下的贝叶斯滤波问题。为了处理非线性和非高斯滤波问题，20世纪90年代人们提出了粒子滤波19-22,71。粒子滤波给出了非线性和非高斯情况下最优贝叶斯滤波的近似数值解，当前已成为最优贝叶斯滤波理论领域的研究热点。2.2 卡尔曼系滤波器为了实现系统状态估计，在20世纪40年代，Wiener和Kolmogorov彼此独立地提出了一种最优线性滤波方法，称为维纳
18、滤波78。但维纳滤波的不足之处是：（1）求解比较复杂，很难得到解析解，而且求得的最优滤波器在工程上很难实现；（2）计算是非递推的，需存储全部观测数据，存储量大且实时性差；（3）不适用于非平稳随机过程的滤波。在20世纪60年代，Kalman突破了维纳滤波的局限性，提出了在时域上的状态空间方法。引入状态变量和状态空间概念，从而改变了对滤波问题的一般描述，即它不是要求直接给出信号过程的二阶特性或谱密度函数，而是把信号过程视为在白噪声作用下一个线性系统的输出。在此基础上，利用Hilbert空间投影理论，Kalman提出了状态估计理论，称为卡尔曼滤波理论74。卡尔曼滤波器给出了一套易于在计算机上实时实现
19、的最优线性递推滤波算法，适合处理多变量系统和时变系统，适合处理非平稳随机过程，克服了维纳滤波理论的缺点和局限性70,78。但是，卡尔曼滤波器只能有效处理线性高斯情况下的滤波问题。在滤波器理论中，扩展卡尔曼滤波是应用最为广泛的一种非线性状态估计方法，其实质是将非线性的状态转移模型或观测模型线性化，然后利用卡尔曼预测方程和更新方程递推求解79,80。然而，非线性系统状态空间模型的线性化经常会在状态估计中引入较大误差。为了更好地处理非线性状态估计问题，Julier等在扩展卡尔曼滤波的理论框架下提出了Unscented卡尔曼滤波52,53。在此，我们将卡尔曼滤波、扩展卡尔曼滤波和Unscented卡尔
20、曼滤波统称为卡尔曼系滤波器。本节将在最优贝叶斯滤波理论框架下简要介绍卡尔曼系滤波器理论。2.2.1 经典卡尔曼滤波器设动态系统S具有线性、高斯性，则式（2.14）定义的状态空间模型可改写为： (2.23)其中，为系统状态转移矩阵，由线性状态转移模型确定；为系统观测矩阵，由线性观测模型确定；系统噪声和观测噪声为零均值高斯白噪声（即和服从高斯分布），且其协方差阵为：，。令系统状态初始分布为高斯分布，可以证明，若服从高斯分布 (2.24)则和也是高斯的， (2.25) (2.26)其中，为高斯分布函数，为均值，为协方差。于是，通过式（2.24）-（2.26）的卡尔曼滤波递推求解，即可实现式（2.21
21、）和式（2.22）所定义的贝叶斯滤波递推求解70。卡尔曼滤波可分为两部分：一步预测和观测更新。其算法如下：算法 2.1（卡尔曼滤波算法）：（1）初始化：对于，给定初始状态的均值和方差；（2）递推求解：对于，则有（a）一步预测： (2.27) (2.28)（b）观测更新 (2.29) (2.30) (2.31)由于式（2.24）-（2.26）可由其均值和协方差阵完全表征，因此卡尔曼滤波是递推贝叶斯最优滤波的显式解。也就是说，在给定线性高斯假设条件下，卡尔曼滤波与贝叶斯最优滤波是完全等价的。2.2.2 扩展卡尔曼滤波器对于式（2.14）定义的系统状态空间模型，若状态转移模型和状态观测模型是非线性的
22、（但仍假设具有零均值高斯白噪声），则经典卡尔曼滤波不能实现最优滤波问题。一般地，人们通过各种非线性近似求得近似解。最基本的近似方法是泰勒近似法，其思路是：当状态的先验分布可用高斯分布近似时，状态的条件分布完全由其均值和协方差阵表征，若在状态的滤波值和预测值周围分别将状态转移方程和观测方程进行泰勒展开： (2.32) (2.33)其中，和为被截掉的二阶以上的高阶项。此种近似称为局部线性近似（也称为一阶泰勒近似），和是线性化的雅可比阵。 (2.34)于是，非线性高斯系统S的局部线性化系统为 (2.35)其中，在给定和时，和分别为确定性分量；和为随机噪声，且包含了线性化误差，因而已是非高斯的。如果忽
23、略线性化误差，即和，那么式（2.35）所示的局部线性化系统具有线性、高斯模型，从而可利用2.2.1节的线性高斯系统的卡尔曼滤波公式递推求解，即可近似假定状态的条件分布是高斯的： (2.36) (2.37) (2.38)显然，式（2.36）-（2.38）中的参数可通过卡尔曼滤波公式递推近似求解。此种非线性滤波方法称为扩展卡尔曼滤波（EKF），其算法如下：算法 2.2（EKF算法）：（1）初始化：对于，给定状态先验高斯分布的均值和方差；（2）递推求解：对于，则有（a）一步预测： （b）观测更新 从EKF的机理分析可知，EKF是一种局部次优的贝叶斯滤波估计，且当系统的非线性较强、状态的条件分布用高斯
24、分布近似的误差较大时，采用EKF近似非线性滤波可能导致较大的累积估计误差。一般地，EKF在应用中要注意两点69,70,77,78：（1）基于泰勒展开的线性化方法易受参考点的影响。EKF是在当前估计值处进行泰勒展开的，并取其线性近似。在EKF递推计算过程中，卡尔曼滤波增益依赖于当前的状态估计值。如果当前估计值与真实值相差很大，则参考点的偏离将引起进一步的线性化误差以及不精确的卡尔曼滤波更新。（2）由于EKF使用了两个雅可比矩阵的计算，所以在EKF应用时应注意状态转移模型和观测模型的连续性。以上两点构成了EKF的基本应用前提：小偏差初始条件和系统较弱的非线性；且和足够光滑，以确保雅可比阵、的存在性
25、。2.2.3 Unscented卡尔曼滤波器EKF是在滤波器理论中应用较为广泛的一种非线性状态估计方法，但当非线性程度比较高时，非线性模型的线性化将导致较大的误差。为了有效地处理非线性状态估计问题，Julier等提出了Unscented卡尔曼滤波（UKF）52,53,81。UKF的核心是Unscented变换，其使用一组适当选择的加权的离散采样（又称为Sigma点）表征系统状态概率分布的均值和协方差，这组采样点根据非线性系统状态空间模型进行预测和观测，从而不需线性化81,82。1. Unscented变换与尺度Unscented变换Unscented变换是一种近似计算经历了非线性变换的随机变量
26、的统计特性的新方法，它建立的动机是：近似一个概率分布比近似一个任意的非线性变换（或函数）要容易81-83。假设是均值为、协方差阵为的维随机向量，且经历非线性变换： (2.39)其中，为非线性函数。为了估计随机向量的均值和协方差阵，Julier和Uhlmann提出了Unscented变换（UT）52,82。其基本思想是：根据一种特定的、确定性的方法采样个加权随机样本，这组随机样本的均值为、协方差阵为；然后将非线性函数作用于每个样本得到变换后的一组随机样本，且能很好地表征随机变量的统计特性。UT的具体步骤如下：（1）根据下列方程采样个加权随机样本点， (2.40) (2.41)其中，是样本点权，且
27、；是尺度因子，控制采样点与样本均值之间的距离；是的第列（或行）的矩阵平方根。这些加权样本点称为Sigma点，并记为。（2）将非线性函数作用于每个Sigma点，则得到一组随机样本点： （3）估计随机向量的均值和协方差阵， (2.42)显然地，UT不需要将非线性函数线性化，也不需要计算雅可比矩阵。而且可以证明52：UT能精确估计任意非线性函数的二阶泰勒近似解；估计误差为三阶及三阶以上高阶矩项截断误差，且该估计误差被尺度化。在UT中，Sigma点数随着状态空间的维数的增大而增大；且Sigma点在状态空间的分布情况决定了UT的性能。特别地，在严重的非线性情况下Sigma点的在状态空间的分布是影响UT性
28、能的关键。针对这个问题，Sigma点被尺度化分布在状态空间里，第个Sigma点到均值的距离为，且距离尺度化比为。但是，当时，权使得估计的协方差阵可能是非半正定的。鉴于此，Julier提出了尺度化Unscented变换（SUT），其不仅能保持估计的二阶精度，而且能使协方差估计是半正定的82。在相同计算代价条件下，SUT能部分地引入高阶矩信息，从而提高估计的精度。对于SUT，Sigma点按如下策略进行采样： (2.43) (2.44) (2.45)其中，是均值估计权，是协方差阵估计权；是尺度化因子，将控制Sigma点的分布；是加权项，将引入高阶矩信息，提高估计精度（高斯先验下最优值为）。于是，随机
29、向量的均值和协方差阵可按下式计算： (2.46)2. Unscented卡尔曼滤波器Julier等提出将SUT和卡尔曼滤波结合实现高斯、非线性情况下的贝叶斯递推滤波问题，这种滤波称为Unscented卡尔曼滤波器（UKF）53,81,83。如果假设式（2.14）定义的系统状态空间模型是高斯、非线性的，并将噪声变量引入到状态变量中产生扩展的状态变量，则UKF算法具体如下：算法 2.3（UKF算法）：（1）初始化：对于，令， (2.47)其中，和是初始状态的均值和协方差阵；系统噪声和观测噪声为零均值高斯白噪声，且其协方差阵为和。（2）递推估计：对于，则有（a）计算Sigma点 (2.48)（b）U
30、KF预测 (2.49) (2.50) (2.51) (2.52) (2.53)（c）UKF更新 (2.54) (2.55) (2.56) (2.57) (2.58)其中，是扩展状态空间的维数，是卡尔曼增益矩阵。与EKF相比，UKF实现高斯、非线性滤波不需计算雅可比阵，能实现任意高斯、非线性情况下的状态估计问题。无论是在理论上还是在实际应用中，UKF被证明都要优于EKF52,53,81,83。2.3 粒子滤波对于式（2.21）和式（2.22）定义的递推贝叶斯滤波问题，卡尔曼系滤波器在对状态后验分布递推求解中，仅对后验分布的一阶和二阶矩进行近似递推计算。因此，卡尔曼系滤波器仅适用于高斯、线性系统或
31、高斯、非线性系统的滤波问题。对于非高斯、非线性系统，其状态分布实际上都有无穷个参数，于是仅在递推参数中传递两个低阶矩参数（或高斯分布假设）是不够的。显然地，非参数估计方法是一种有效的途径，将可完全放弃对状态分布所作的高斯假设。粒子滤波就是这样一种方法，采用序贯蒙特卡罗（Sequential Monte Carlo）模拟来近似状态分布，并实现贝叶斯递推滤波。粒子滤波器（又称为CONDENSATION、Bootstrap Filter或Sequential Monte Carlo Filter）分别由信号处理21,84、计算机视觉20,85、统计22,71等领域独立地提出用以解决非高斯、非线性贝叶
32、斯递推滤波问题。粒子滤波是以Monte Carlo随机模拟理论为基础，将系统状态后验分布用一组加权随机样本（称为粒子）近似表示，新的状态分布通过这些随机样本的贝叶斯递推估计。粒子滤波主要包括三个步骤：（1）粒子采样，从建议分布（Proposal Distribution）中抽取一组新的粒子；（2）粒子加权，根据观测概率分布和贝叶斯公式计算各个粒子的权值；（3）估计输出，输出系统状态的均值、协方差或高阶矩等。此外，为了避免粒子滤波中出现的退化现象，重采样步骤经常被采用。本节首先简单介绍蒙特卡罗随机模拟原理，在此基础上阐述标准粒子滤波的理论框架；并讨论标准粒子滤波的粒子退化和“样贫”问题；同时，简
33、单讨论了粒子滤波的收敛性问题。2.3.1 蒙特卡罗（Monte Carlo）随机模拟在很多复杂的统计问题（比如在信号处理、计算机视觉领域的高维贝叶斯推理和组合计算等）中，很难直接进行理论分析并进而求解，而随机模拟是非常实用的方法。随机模拟就是设法按问题的要求与条件去构造一系列的随机样本，并用这些样本的频率代替对应的概率作统计分析和推断86,87。在概率论发展初期，随机模拟原型常常来自博彩，于是人们就以博彩之都蒙特卡罗（Monte Carlo）作为随机模拟别称。蒙特卡罗方法经常应用于求解高维积分和优化问题，广泛应用于物理学88、统计学22,71,89-91、信号处理21,25,84、机器学习87
34、,92以及计算机视觉20,93-96等领域。下面简要介绍蒙特卡罗方法的基本原理86,87：对于高维空间上高维积分， (2.59)其中，为定义在高维空间上的概率分布；是关于任意可积函数（多数为非线性函数），且满足。如果从概率分布独立地抽取个随机样本，则样本集是独立同分布的，于是概率分布即可由这些样本近似： (2.60)其中，为在样本处的delta-Dirac函数。于是，式（2.59）定义的高维积分问题可近似为如下求和问题： (2.61)这种基于随机模拟的积分方法称为蒙特卡罗积分。由式（2.61）可知，蒙特卡罗积分是几乎处处收敛于；而且，由大数定律可知，是无偏的。如果的方差是有界的，且定义，则的方
35、差为，并且由中心极限定理则有： (2.62)如果概率分布具有标准形式（比如高斯分布），那么可直接从采样获得随机样本集87。但是，在实际应用中很难找到标准形式，于是很多复杂的采样策略被提出，比如拒绝采样法97,98、重要性采样法49,60,86,97,99和马尔科夫链-蒙特卡罗采样法86,87,94,95,97等。在2.4节中将详细讨论随机模拟的各种采样策略问题。2.3.2 标准粒子滤波器对于式（2.20）定义的递推贝叶斯状态估计问题，系统的后验分布可以由蒙特卡罗方法近似。在非高斯、非线性情况下，后验分布很难具有标准形式，因此不能直接从后验分布采样。如果令为定义在非高斯、非线性系统状态空间上易于
36、抽样的条件概率分布，且与后验分布具有相同或更大的支撑集 函数在其定义域全体函数值不为零的点的集合称为的支撑集，记为。，则称为建议分布（也称为重要性函数）49,86,87。设从建议分布抽取个随机样本形成样本集，则有： (2.63)于是，对于后验分布则有： (2.64)其中，为归一化的样本权，且，其中 (2.65)由此可见，后验概率可由一组加权的随机样本（称为粒子）近似，而这样的采样方法，称为重要性采样方法49,86,87。为了实现递推贝叶斯状态估计，可将建议分布改写成递推形式。如果假设当前系统状态独立于未来观测值，则有： (2.66)并且，由递推贝叶斯状态估计公式（式（2.20）可知： (2.6
37、7)如果将式（2.66）和式（2.67）代入式（2.65），对于粒子，其权值即为： (2.68)为了方便推导，式（2.67）未引入归一化因子。显然地，在式（2.67）引入归一化因子等价于将权值归一化。如果假设状态是一马尔科夫过程且和观测是条件独立的，那么对于递推贝叶斯滤波问题，粒子的权可根据式（2.68）改写为： (2.69)于是，后验概率可近似为： (2.70)特别地，如果建议分布为状态先验分布，并递推化可得： (2.71)将式（2.71）代入式（2.69），则有： (2.72)一般地，把这样的递推贝叶斯滤波称为标准粒子滤波41,49,61。标准粒子滤波算法如下：算法 2.4（标准粒子滤波算
38、法）：（1）初始化：对于，根据状态先验分布建立初始状态粒子集，其中；（2）对于（a）采样：从状态转移概率采样得到新的粒子集；（b）权值计算：根据式（2.72）计算粒子的权，并进行归一化；（3）状态估计输出：根据粒子集，状态后验分布可由式（2.70）近似，且状态估计为；（4）重采样：根据粒子的权从粒子集重新抽取个粒子，并令，则建立新的粒子集。2.3.3 退化现象和“样贫”问题在标准粒子滤波中，经常出现退化现象，其表现为：经过若干次递推计算后，除了少数粒子外，其余粒子的权值可忽略不计，从而使得大量递推计算浪费在对几乎不起任何作用的粒子的更新上，甚至最后只剩下一个权值很大（几乎接近1）的有效粒子，而
39、其他粒子的权值几乎为零，从而产生一个退化分布48,49。粒子的退化现象产生于有限抽样和计算机的有限字长。当粒子的建议分布与后验分布具有窄的支撑集，且重叠部分很少时，由于有限抽样，大部分粒子不在公共支撑集上，因此经测量的贝叶斯更新后，将具有较小的权值。经过多次递推计算以及计算机下溢，导致大部分粒子的权值为零，从而成为无效粒子。为了度量粒子集的退化，Liu等97,100提出了有效粒子数，且定义为： (2.73)由上式可知，很难精确计算，其近似值可通过下式计算： (2.74)显然地，一个或几个有效粒子是不能近似状态的统计特性（或概率分布）的，因此退化现象将严重影响粒子滤波的性能。一般地，避免退化现象
40、的方法有两种：一种方法是选取有效的建议分布法（在2.4节讨论）；另一种方法就是重采样法。重采样法就是在粒子滤波算法的权值计算步骤之后，引入重采样步骤（见算法2.4）。重采样就是根据粒子的权值对粒子集进行重新采样，权值大的粒子被多次抽取（抽样次数记为，且），而权值小的粒子被随机剔除48,49。在粒子滤波中，常用的重采样算法为多项采样算法71。对于多项采样算法，其采样方差为。为了减小重采样的方差和提高采样的有效性，残差采样71,101和系统采样22,102等重采样算法被提出，而系统采样是一种采样方差较小的重采样算法102。基于系统采样的重采样算法如下：算法 2.5（多项重采样算法）：对于粒子集，则
41、有（1）令，并计算累积权值；（2）产生上均匀分布的随机数：；（3）令，对于，（a）；（b）如果，则；否则，；虽然重采样可以消除小权值粒子在粒子滤波中的影响，但是又引入了新的负面问题，称为“样贫”49。特别地，当系统噪声较小时，“样贫”问题将尤为严重。“样贫”问题源于重采样过程。由于重采样将高权值的粒子过度复制，有效粒子数在重采样后将减少，导致新的粒子集的信息容量降低。这样一来，经过若干次递推计算以后，有效粒子几乎被耗尽，粒子集的信息容量将严重降低，很难反映真实的状态统计特性。减轻“样贫”问题的常用方法有：重采样移动法103、正则化法104和马尔科夫链-蒙特卡罗法105,106等。2.3.4粒子滤波的收敛性问题粒子滤波是一种有效的非高斯、非线性递推贝叶斯滤波方法，而且粒子滤波又是一种全