粒子滤波器

✍ dations ◷ 2024-12-23 04:24:34 #蒙地卡罗方法,计算统计学,控制理论,机器人控制,统计力学

粒子滤波器(particle filter)是一种使用蒙特卡罗方法(Monte Carlo method)的递归滤波器，透过一组具有权重的随机样本(称为粒子)来表示随机事件的后验概率，从含有噪声或不完整的观测序列，估计出动力系统的状态，粒子滤波器可以运用在任何状态空间的模型上。

粒子滤波器是卡尔曼滤波器(Kalman filter)的一般化方法，卡尔曼滤波器建立在线性的状态空间和高斯分布的噪声上；而粒子滤波器的状态空间模型可以是非线性，且噪声分布可以是任何型式。

从统计和概率的角度来看，粒子滤波器属于分支/遗传算法的类别，并且是指平均场类型相互作用的粒子方法。这些粒子方法的解释取决于科学专业。在进化计算中，平均场遗传类型粒子方法通常用作启发式和自然搜索算法（也称为元启发算法）。在计算物理学和分子化学中，它们用于解决费曼-卡兹(Feynman-Kac)路径积分问题，或者用于计算玻尔兹曼-吉布斯(Boltzmann-Gibbs)度量，薛定谔(Schrödinger)算子的最高特征值和基态。在生物学与遗传学中，它们还代表了某些环境中一群个体或基因的进化。

平均场类型进化计算技术的起源可以追溯到1950年和1954年，这是艾伦·图灵(Alan Turing)在遗传类型突变选择学习机上的开创性工作，以及新泽西州普林斯顿高级研究所的尼尔斯·艾尔·巴里切利(Nils Aall Barricelli)的文章。统计学中的粒子滤波器最早可追溯到1950年代中期。哈默斯利(Hammersley)等人在1954年提出的“穷人的蒙特卡洛”法暗示了当今使用的遗传类型粒子过滤方法。 1963年，尼尔斯·艾尔·巴里切利(Nils Aall Barricelli)模拟了一种遗传类型算法，以模仿个人玩简单游戏的能力。在进化计算文献中，遗传类型突变选择算法通过1970年代初期约翰·霍兰德(John Holland)的开创性工作而流行，特别是1975年出版的著作。

在生物学与遗传学中，澳大利亚遗传学家亚历克斯·弗雷泽（Alex Fraser）在1957年还发表了一系列有关人工选择生物的遗传类型模拟的论文。由生物学家进行的进化计算机模拟在1960年代初变得更加普遍，并且该方法在弗雷泽(Fraser)和伯内尔(Burnell)（1970）和克罗斯比(Crosby)（1973）的书中进行了描述。弗雷泽(Fraser)的模拟包括了现代突变选择遗传粒子算法的所有基本要素。

从数学观点来看，给定一些部分和噪声观测结果的信号随机状态的条件分布，是由信号的随机轨迹上的费曼-卡兹(Feynman-Kac)概率描述的，该概率由一系列似然势函数加权。量子蒙特卡洛法，更具体地说是扩散蒙特卡洛法，也可以解释为费曼-卡兹路径积分的平均场遗传类型粒子近似。量子蒙特卡罗方法的起源通常归因于恩里科·费米(Enrico Fermi)和罗伯特·里克特迈耶(Robert Richtmyer)，他们于1948年开发了一种中子链反应的均值场粒子解释方法，但它是第一个类似启发式和遗传类型的粒子算法（又名“重采样”或“重新配置”蒙特卡洛方法））是在1984年由Jack H. Hetherington提出的用于估计量子系统（在简化的矩阵模型中）的基态能量的方法。还可以引用1951年出版的西奥多·爱德华·哈里斯(Theodore E. Harris)和赫尔曼·康恩(Herman Kahn)的开创性著作，该著作使用平均场，但类似于启发式遗传方法，用于估计粒子传输能量。在分子化学中，遗传启发式粒子方法（又名修剪和富集策略）的使用可以追溯到1955年马歇尔（Marshall）的开创性工作。罗森布鲁斯（N. Rosenbluth）和阿里安娜·W·罗森布卢斯(Arianna. W. Rosenbluth)。

遗传粒子算法在高级信号处理和贝叶斯推理中的使用是最近的。 1993年1月，北川源四郎开发了“蒙特卡罗滤波器”，该文在1996年出现了轻微修改。1993年4月，戈登等人在他们的开创性著作中发表了遗传类型算法在贝叶斯统计推断中的应用。作者将其算法命名为“自举滤波器”，并证明与其他滤波方法相比，其自举算法不需要对该状态空间或系统噪声进行任何假设。独立的是皮埃尔·德尔·莫拉尔(Pierre Del Moral)和希米尔孔·卡瓦略(Himilcon Carvalho)，皮埃尔·德尔·莫拉尔(Pierre Del Moral)，安德烈·莫宁(André Monin)和热拉尔·萨鲁特(Gérard Salut)在1990年代中期发表的有关粒子过滤器的论文。 1989-1992年初，皮埃尔·德尔·莫拉尔(P.Del Moral)，J.C.诺耶(J.C. Noyer)，G.Rigal和G.Salut在LAAS-CNRS中也开发了信号处理技术，并通过STCAN（服务技术）进行了一系列受限和分类研究des Constructions and Armes Navales），IT公司DIGILOG和LAAS-CNRS（系统分析和体系结构实验室）来解决RADAR / SONAR和GPS信号处理问题。

从1950年到1996年，所有有关粒子过滤器，遗传算法的出版物，包括在计算物理和分子化学中引入的修剪和重采样蒙特卡洛方法，都提出了适用于不同情况的类似于自然和启发式算法，而没有任何证据证明它们的一致性，也没有讨论估计的偏差以及基于族谱和祖先树的算法。

这些粒子算法的数学基础和首次严格的分析是由皮埃尔·德尔·莫拉尔(Pierre Del Moral)在1996年提出的。本文还包含了似然函数和未标准化的条件概率测度的粒子近似的无偏性质的证明。本文介绍的似然函数的无偏粒子估计器如今已用于贝式统计推理中。

到1990年代末，Dan Crisan，Jessica Gaines和Terry Lyons以及Dan Crisan，Pierre Del Moral和Terry Lyons也开发了具有不同种群大小的分支类型粒子方法。 P. Del Moral，A。Guionnet和L. Miclo在2000年开发了该领域的进一步发展。第一个中心极限定理是由Pierre Del Moral和Alice Guionnet于1999年以及Pierre Del Moral和Laurent Miclo于2000年提出的。Pierre于1990年代末提出了关于粒子过滤器时间参数的第一个统一收敛结果。 Del Moral和Alice Guionnet。 2001年，P。Del Moral和L. Miclo首次对基于族谱树的粒子滤波器平滑器进行了严格的分析。

关于Feynman-Kac粒子方法论的理论和相关的粒子过滤器算法已在2000年和2004年的书中得到发展。这些抽象的概率模型封装了遗传类型算法，粒子和自举滤波器，交互的卡尔曼滤波器（又称Rao–Blackwellized粒子滤波器），重要性采样和重采样样式的粒子滤波器技术，包括基于族谱树的方法和用于解决滤波和平滑问题的粒子后向方法。其他类别的粒子滤波方法包括基于族谱树的模型，后向马尔可夫粒子模型，自适应平均场粒子模型，岛型粒子模型和粒子马尔可夫链蒙特卡罗方法。

粒子滤波器的目标是通过观测值估计状态变量的后验概率。粒子滤波器是针对隐马尔可夫模型设计的，系统中既包含隐藏的变量，也包含可观测的变量。可观测变量（观测过程）与隐藏变量（状态过程）通过某种已知的函数关系相关联。类似地，描述状态变量演化的动态系统也是概率学上已知的。

一个通用粒子滤波器通过观察测量过程估计隐藏状态的后验分布。考察下图所示的状态空间：

滤波问题是要通过给定的观测过程 $Y_{0},\cdots ,Y_{k},$ 的值，序列地估计隐藏状态 $X_{k}$ ( | ₀,₁,…,) 。粒子滤波器方法使用经验测量和通用粒子算法，提供了这些条件概率的近似值。与之对比，马尔可夫链蒙特卡洛（MCMC）或重要性采样方式会对整个后验概率 (₀,₁,…, | ₀,₁,…,) 进行建模。

粒子滤波器能够从一系列含有噪声或不完整的观测值中，估计出动力系统的内部状态。在动力系统的分析中，需要两个模型，一个用来描述状态随时间的变化(系统模型)，另一个用来描述每个状态下观测到的噪声(观测模型)，将这两个模型都用概率来表示。

在许多情况下，每得到一个新的观测值时，都必须对系统做出一次估计，利用递归滤波器，能够有效地达到这样的目的。递归滤波器会对得到的资料做连续处理，而非分批处理，因此不需要将完整的资料储存起来，也不需要在得到新的观测值时，将现有的资料重新做处理。递归滤波器包含两个步骤：

预测：利用系统模型，由前一个状态的资讯预测下一个状态的概率密度函数。

更新：利用最新的观测值，修改预测出的概率密度函数。

借由贝叶斯推论(Baysian inference)，我们可以描述出状态空间的概率，并在得到新的观测值时，对系统做出更新，因而达成上述目的。

在某些问题中，在信号的随机状态下，观测值的条件分布可能不具有密度，或者不可能或太复杂而无法计算。在这种状况下，我们需要求近似值。其中一个策略是借由马尔可夫链 ${\mathcal {X}}_{k}=\left(X_{k},Y_{k}\right)$ ${\mathcal {X}}_{k}=\left(X_{k},Y_{k}\right)$ 来取代信号 $X_{k}$ $X_k$ ，采用虚拟观测的形式

对于具有已知概率密度函数的独立序列中的某些序列，核心的概念是观测此式子(公式)

在部分的观测 ${\mathcal {Y}}_{0}=y_{0},\cdots ,{\mathcal {Y}}_{k}=y_{k},$ ${\mathcal {Y}}_{0}=y_{0},\cdots ,{\mathcal {Y}}_{k}=y_{k},$ 之下，粒子滤波器是和马尔可夫链 ${\mathcal {X}}_{k}=\left(X_{k},Y_{k}\right)$ ${\mathcal {X}}_{k}=\left(X_{k},Y_{k}\right)$ 有相关的，定义为根据演化出的 $\mathbb {R} ^{d_{x}+d_{y}}$ $\mathbb {R} ^{d_{x}+d_{y}}$ ，在一些明显的滥用性符号 $p({\mathcal {Y}}_{k}|{\mathcal {X}}_{k})$ $p({\mathcal {Y}}_{k}|{\mathcal {X}}_{k})$ 下，带有似然函数的粒子。这些概率性的技术跟近似贝式计算非常相关。在量子滤波器中，这些近似贝式量子滤波器的技术在1998年被皮埃尔·德尔·莫拉尔(P. Del Moral), 让·雅各德(J. Jacod)和菲力普·普罗特(P. Protter)所采用。且被皮埃尔·德尔·莫拉尔, 阿诺·杜斯(A. Doucet)和阿杰·贾斯拉(A. Jasra)更进一步发展。

借由贝式定理在条件概率中可以得出

当

粒子滤波器也是一个近似值，当有足够的粒子时，结果会更加的准确。非线性滤波方程可以借由递归得出

${\begin{aligned}p(x_{k}|y_{0},\cdots ,y_{k-1})&{\stackrel {\text{updating}}{\longrightarrow }}p(x_{k}|y_{0},\cdots ,y_{k})={\frac {p(y_{k}|x_{k})p(x_{k}|y_{0},\cdots ,y_{k-1})}{\int p(y_{k}|x'_{k})p(x'_{k}|y_{0},\cdots ,y_{k-1})dx'_{k}}}\\&{\stackrel {\text{prediction}}{\longrightarrow }}p(x_{k+1}|y_{0},\cdots ,y_{k})=\int p(x_{k+1}|x_{k})p(x_{k}|y_{0},\cdots ,y_{k})dx_{k}\end{aligned}}$ ${\begin{aligned}p(x_{k}|y_{0},\cdots ,y_{k-1})&{\stackrel {\text{updating}}{\longrightarrow }}p(x_{k}|y_{0},\cdots ,y_{k})={\frac {p(y_{k}|x_{k})p(x_{k}|y_{0},\cdots ,y_{k-1})}{\int p(y_{k}|x'_{k})p(x'_{k}|y_{0},\cdots ,y_{k-1})dx'_{k}}}\\&{\stackrel {\text{prediction}}{\longrightarrow }}p(x_{k+1}|y_{0},\cdots ,y_{k})=\int p(x_{k+1}|x_{k})p(x_{k}|y_{0},\cdots ,y_{k})dx_{k}\end{aligned}}$

(Eq. 1)

依照惯例 $p(x_{0}|y_{0},\cdots ,y_{k-1})=p(x_{0})$ $p(x_{0}|y_{0},\cdots ,y_{k-1})=p(x_{0})$ 在k=0时。非线性滤波方程问题的关键在于依序计算这些条件概率分布。

固定一个时间范围和一个序列的观测 $Y_{0}=y_{0},\cdots ,Y_{n}=y_{n}$ $Y_{0}=y_{0},\cdots ,Y_{n}=y_{n}$ ，在k=0....n，我们定义:

在这种表示法中，对于从原点k = 0到时间k = n的轨迹集合上的任何有界函数F，我们都可以用费曼-卡兹公式来表示

这些费曼-卡兹路径积分模型出现在各种科学学科中，包括计算物理，生物学，信息论和计算机科学。他们的解释取决于应用领域。举例来说，如果我们选择状态空间某些子集的指标函数 $G_{n}(x_{n})=1_{A}(x_{n})$ $G_{n}(x_{n})=1_{A}(x_{n})$ ，它们代表马尔可夫链在给定管中的条件分布。我们会得到:

和

只要标准化常数严格为正。

在动力系统中，我们常常需要对物体做追踪。假设有一动力系统，已知其状态序列 $\{x_{k},k\in N\}$ $\{x_{k},k\in N\}$ 定义为

$x_{k}=f_{k}(x_{k-1},v_{k-1})$ $x_{k}=f_{k}(x_{k-1},v_{k-1})$

其中 $f_{k}:R^{n_{x}}\times R^{n_{v}}\rightarrow R^{n_{x}}$ $f_{k}:R^{n_{x}}\times R^{n_{v}}\rightarrow R^{n_{x}}$ 是状态转移函数，可以是非线性的函数， $\{v_{k-1},k\in N\}$ $\{v_{k-1},k\in N\}$ 是一个独立且相同分布(i.i.d.)的过程噪声序列， $n_{x}$ $n_x$ 和 $n_{v}$ $n_{v}$ 分别代表状态和过程噪声向量的维度， $N {\displaystyle N}$ $N$ 为自然数的集合。追踪的目标是要递归地从观测值 $y_{k}$ $y_{k}$ 估计出 $x_{k}$ $x_k$ ，而观测值 $y_{k}$ $y_{k}$ 定义为

$y_{k}=h_{k}(x_{k},n_{k})$ $y_{k}=h_{k}(x_{k},n_{k})$

其中 $h_{k}:R^{n_{x}}\times R^{n_{n}}\rightarrow R^{n_{y}}$ $h_{k}:R^{n_{x}}\times R^{n_{n}}\rightarrow R^{n_{y}}$ 是观测函数，可以是非线性的函数， $\{n_{k},k\in N\}$ $\{n_{k},k\in N\}$ 是一个独立且相同分布(i.i.d.)的观测噪声序列， $n_{y}$ $n_y$ 和 $n_{n}$ $n_{n}$ 分别代表观测值和观测噪声向量的维度。

从贝叶斯理论的观点来看，追踪问题是要根据到时间 $k {\displaystyle k}$ $k$ 为止的已知资讯 $y_{1:k}$ $y_{1:k}$ ，递归地计算出时间 $k {\displaystyle k}$ $k$ 的状态 $x_{k}$ $x_k$ 的可信度，这个可信度用概率密度函数 $p(x_{k}\mid y_{1:k})$ $p(x_{k}\mid y_{1:k})$ 来表示。假设初始概率密度函数 $p(x_{0}\mid y_{0})\equiv p(x_{0})$ $p(x_{0}\mid y_{0})\equiv p(x_{0})$ (称为先验概率)为已知，则利用“预测”和“更新”两个步骤就能递归地计算出概率密度函数 $p(x_{k}\mid y_{1:k})$ $p(x_{k}\mid y_{1:k})$ 。

在此假设在时间 $k-1$ $k-1$ 的概率密度函数 $p(x_{k-1}\mid y_{1:k-1})$ $p(x_{k-1}\mid y_{1:k-1})$ 为已知。

利用查普曼-科尔莫戈罗夫等式(Chapman–Kolmogorov equation)，可以由状态转移函数与时间 $k-1$ $k-1$ 的概率密度函数 $p(x_{k-1}\mid y_{1:k-1})$ $p(x_{k-1}\mid y_{1:k-1})$ ，计算出时间 $k {\displaystyle k}$ $k$ 的先验概率 $p(x_{k}\mid y_{1:k-1})$ $p(x_{k}\mid y_{1:k-1})$

${\begin{aligned}p(x_{k}\mid y_{1:k-1})&=\int p(x_{k},x_{k-1}\mid y_{1:k-1})d{x_{k-1}}\\&=\int p(x_{k}\mid x_{k-1},y_{1:k-1})p(x_{k-1}\mid y_{1:k-1})d{x_{k-1}}\\&=\int p(x_{k}\mid x_{k-1})p(x_{k-1}\mid y_{1:k-1})d{x_{k-1}}\\\end{aligned}}$