Sqare Root SAM

Simultaneous Localization and Mapping via Square Root Information Smoothing 基于平方根信息平滑法的同步定位和映射

后面的其他部分 和 因子图优化相关问题重了

摘要

我们研究平滑方法作为一个可行的替代扩展卡尔曼滤波的解决方案。特别地,我们研究了将相关信息矩阵或度量雅可比矩阵分解为平方根形式的方法。与EKF相比,这种技术有几个显著的优势:它们更快又精确,可以批量或增量模式使用,可以更好地处理非线性过程和测量模型,并以更低的成本产生整个机器人轨迹。 此外,以一种间接但引人注目的方式,列排序启发式自动利用了SLAM问题的地理本质中固有的局部性。在本文中,我们提出了这些方法背后的理论,以及根据与SLAM问题相关的图形模型对因子分解的解释。我们在大规模环境中展示了模拟结果和实际的SLAM实验,这强调了这些方法作为基于EKF的方法的替代方案的潜力。

1.介绍

其中最早和最流行的方法是基于扩展的卡尔曼滤波器(EKF)[68,60,59,2,69,51,49]。EKF递归地估计机器人当前姿态和所有地标的位置(地图)上的高斯密度。然而,众所周知,EKF的计算复杂度很快就变得相当棘手,因此大量的工作都集中在修改和扩展过滤方法,以应对更大规模的环境。[62, 17, 21, 8, 44, 47,37, 63, 48, 5, 75, 38, 65].然而,当应用于固有的非线性SLAM问题[43]时,过滤本身已经被证明是不一致的,即,即使是占用大量实验的平均值也偏离了真正的解。由于这主要是由于在过滤框架中无法撤销的线性化选择,最近人们对SLAM问题的平滑版本产生了相当大的兴趣

Smoothing and Mapping or SAM

SLAM的平滑方法不仅涉及最新的机器人位置,还涉及整个机器人的轨迹。许多作者考虑了平滑机器人轨迹的问题,这特别适用于传感器,如激光距离查找器,容易在附近的机器人姿态之间产生成对约束。更一般地说,我们可以考虑完整的SLAM问题[74],即最优估计整个传感器姿态集以及环境中所有特征的参数的问题。事实上,这个问题在测量[35]、摄影测量[7,36,66,10]、被称为“束调整”和计算机视觉[24,70,71,76,40]、被称为“运动结构”方面有着悠久的历史。特别是在过去的五年里,有一系列的工作,这些想法被应用到SLAM的背景下。

在本文中,我们展示了平滑是如何成为基于过滤的方法的一种非常快速的替代方法,并且在许多情况下,保持轨迹周围有助于而不是损害性能。特别地,与全SLAM相关的优化问题可以用稀疏线性代数简明地表示,它传统上涉及最小二乘问题[34]的解。在这个框架中,我们研究了将信息矩阵I或度量雅可比矩阵A分解为平方根形式,并应用于同时平滑和映射(SAM)的问题。因为它们是基于矩阵平方根的,所以我们将这类方法称为平方根SAM,或简称√SAM,首先在[18]中引入。我们认为√SAM是SLAM比EKF更好的解决SLAM问题的方法:

  1. 与扩展的卡尔曼滤波协方差或信息矩阵相比,两者都随着时间的推移而完全密集。与平滑相关联的信息矩阵I是和保持稀疏的。
  2. 在典型的映射场景中,这个矩阵I,或者,度量雅可比矩阵A,是映射协方差结构的更紧凑的表示
  3. 两种稀疏,可以分别用稀疏分解或QR分解有效地分解,分别,生成一个平方根信息矩阵R可以用来立即获得机器人的最佳轨迹和地图。

在序列估计文献中,分解信息矩阵被称为平方根信息滤波(SRIF),并于1969年开发用于喷气推进实验室的水手10号金星任务(如[3]所述)。使用协方差或信息矩阵的平方根可以得到更准确和稳定的算法,并且引用梅贝克[57]的话,“许多实践者以相当大的逻辑认为,应该始终采用平方根滤波器,而不是标准的卡尔曼滤波器递归”

梅贝克简要地讨论了SRIF关于平方根过滤的章节,以及它和其他平方根类型的算法是Bierman [3]的一本书的主题。但是,从文献文献数量较少来判断,SRIF和平方根信息平滑器(SRIS)并不常用。

稀疏线性代数,图论,和稀疏√SAM

性能的关键是利用稀疏线性代数中的大量工作,并充分利用与平滑SLAM问题相关的矩阵的稀疏性。性能上最显著的改进来自于在分解矩阵时选择一个良好的变量排序。为了更好地理解这一事实,我们需要研究类似于slam的问题、线性代数和图论之间的密切关系,图论一直是稀疏线性代数30多年的中心[31]。例如,虽然QR分解通常被用作“黑盒”算法,但它实际上是在图上的一种优雅的计算。这是一个反复出现的主题:线性代数的最新进展是数值方法和高级图论的混合,与图形模型文献中常见的推理算法有许多共同的特征。 特别是,稀疏最小二乘解的图论算法是变量消去,其中每个变量(如机器人姿态或地标位置)都用其他变量表示。剔除变量的顺序对QR和Cholesky等矩阵分解算法的运行时间有很大的影响。找到一个最优排序是一个np-完整的问题,但是有一些排序启发式和近似算法在一般问题[1,41]上表现良好,并被内置到像MATLAB这样的程序中。

虽然通用的排序方法极大地提高了性能,但在对变量进行排序时,通过利用SLAM问题的特定图形结构,还可以获得另外15倍的改进。从图表的角度来看SLAM本身有丰富的历史,特别是在过去几年[61,58,29,63,25,30,75,26,74]导致了一些新颖和令人兴奋的发展。下面我们将更仔细地研究图形模型视图和SLAM的稀疏线性代数公式之间的紧密联系。众所周知,信息矩阵I与连接机器人姿态和地标的无向图相关联(见[75])。不太容易理解的事实是,测量雅可比矩阵A是与SLAM相关的因子图的矩阵。此外,平方根信息矩阵R,即分解I或A的结果,本质上对应于一个结树,从图形模型[11]的推理中知道,最近也应用于SLAM [63]。 利用领域知识来获得良好的排序也是线性代数(如[15])的一种趋势,我们相信,通过将这个问题看作是一个图上的计算问题,甚至可以开发出更有效的算法。

2. SLAM and Graph

SLAM是指在定位机器人环境的同时定位机器人的问题,如图1的例子所示。在本节中,我们将介绍SLAM问题,我们使用的符号,并展示文献中已知的三个主要图形模型表示如何对SLAM问题产生独特的视图,从而强调问题的某些方面。在本文的后面,我们建立了这些图和它们的(稀疏)矩阵等价物之间的联系。

下面我们假设我们熟悉基于EKF的SLAM [69,49,8,21]方法。我们不重新推导出扩展的卡尔曼滤波器。相反,在第3节中,我们立即采取了一种平滑的方法,其中地图和机器人的轨迹都被恢复。

SLAM as a Belief Net

信念网也叫贝叶斯网络

根据FastSLAM和其他[61,58,63,29,25,74]设定的趋势,我们通过引用一个信念网络表示来表述这个问题。信念网是一个有向无环图,它编码一组变量的条件独立结构,其中每个变量只直接依赖于图中的前身。我们所采用的模型如图2所示。这里我们用xi表示机器人的时间步长状态,i∈0..M,lj的地标,j∈1..N,用zk测量,用k∈1..K.该网络对应的联合概率模型为

\(P(X,L,Z)=P(x_0)\prod_{i=1}^MP(x_i|x_{i-1},u_i)\prod_{k=1}^KP(z_k|x_{i_k},l_{j_k})\)

其中P(x0)为初始状态的先验,P(xi|xi−1,ui)为运动模型,由控制输入ui参数化,P(z|x,l)为地标测量模型。以上假设在地标l上有一个统一的先验。此外,它假设数据关联问题已经得到解决,即每个测量zk所对应的指标ik和jk都是已知的。

SLAM问题的贝叶斯信念网络表示。机器人的状态x由顶部的马尔可夫链控制,而机器人的环境在底部由一组地标l表示。中间层的测量z由机器人的状态和测量的地标参数控制。 \[ {x_i=f_i(x_{i-1},u_i)+w_i\quad\Leftrightarrow\quad P(x_i|x_{i-1},u_i)\propto\exp-\frac12\|f_i(x_{i-1},u_i)-x_i\|_{\Lambda_i}^2} \] 其中,fi(.)是一个过程模型,以及具有正态分布的具有协方差矩阵Λi的零均值过程噪声,和 \[ z_k=h_k(x_{i_k},l_{j_k})+v_k\quad\Leftrightarrow\quad P(z_k|x_{i_k},l_{j_k})\propto\exp-\frac{1}{2}\|h_k(x_{i_k},l_{j_k})-z_k\|_{\Sigma_k}^2 \] 其中,hk()是一个测量方程,它是正态分布的零平均测量噪声与协方差Σk。在||e||Σ被定义为给定一个协方差矩阵Σ的马氏距离的平方。

上述方程分别模拟了机器人对控制输入及其传感器的响应行为。

这两个式子展示的是状态空间模型中的状态转移的描述,一个是通过随机方程的形式,另一个是通过条件概率密度的形式.

在状态空间模型中,每个状态\(x_i\) 是基于上一个状态\(x_{i-1}\)和一些控制输入\(u_i\),通过一个确定性的函数 \(f_i\) 来预测的,并加上一些随机噪声\(w_i\)。随机噪声通常假定为高斯分布,也就是正态分布,具有零均值和一定的协方差矩阵\(Q_i\)

image-20240510171449639

SLAM as a Factor Gragh

虽然信念网是思考SLAM问题的生成方面的一种非常自然的表示,但因子图与潜在的优化问题有更紧密的联系。由于测量值z_k在图2是已知的(证据,在图形模型的术语中),我们可以自由地消除它们作为变量。相反,我们把它们作为实际未知数上的联合概率因子的参数,这自然导致了众所周知的因子图表示,这是一类二部图图形模型,可以用来表示这种因子密度[46]。在因子图中,有未知数的节点和定义在它们上的概率因子的节点,图的结构表示每个因子涉及哪些未知数。 图1中的示例的因子图如图3所示。可以看出,地标测量zk和测程链接ui都有因子节点.

在讨论SLAM(Simultaneous Localization and Mapping)问题时,我们通常需要处理大量的变量和测量数据。信念网(或贝叶斯网络)和因子图是两种常用的图形表示,用以表达这些变量之间的依赖关系。这两种表示方式各有优缺点,而对于SLAM这类优化问题,因子图提供了一种更直接且与问题求解密切相关的表示形式。

为什么可以自由消除测量值 ( z_k ) 作为变量

在因子图中,我们通常关注的是如何根据已知的测量值(如传感器数据)来推断机器人的位置和地图的配置。测量值 ( z_k ) 通常被视为“已知”的,它们是收集到的传感器数据,不是需要求解的未知数。在构建模型时,我们不需要将这些测量值作为待求解的变量;它们是给定的条件或“证据”,已经存在并用于影响其他未知变量的估计。

为什么我们把它们作为实际未知数上的联合概率因子的参数

在因子图表示中,每一个因子通常表示变量之间的概率关系。例如,在SLAM中,一个因子可能表示在特定的位置和地图配置下获得某个测量值的概率。这样,测量值 ( z_k ) 作为参数出现在因子中,而非作为独立的变量。这使得因子图直接反映了未知变量(如机器人位置和地图配置)之间的条件依赖关系。通过这种方式,因子图将关注点集中在如何根据测量数据来推断或优化未知变量上,而非处理测量数据本身作为未知数。

这自然导致了众所周知的因子图表示

因子图的这种表示方式使得它非常适合于SLAM这类问题的求解,因为SLAM本质上是一个基于大量测量数据对未知状态进行优化的问题。在因子图中,各个因子直接关联到一组变量,并定义这些变量间的局部关系,这对于执行如图搜索、优化和推断等操作是非常有效的。相对于信念网,因子图提供了一种更为紧凑和针对优化问题的自然表示,使得算法能更直接地操作相关的概率分布,从而更有效地实现SLAM。

总之,因子图通过将测量值作为已知的参数而不是变量来简化问题表示,使得模型更直接地关联到优化问题本身,这对于高效解决SLAM等复杂问题至关重要。

因子图是图形模型中的一种,它用于表示多个变量间复杂的依赖关系。在因子图中,节点表示变量(例如,在SLAM中的机器人位置和地图的特征点),而边则表示这些变量之间的条件依赖,即因子。因子通常定义了一组变量之间的概率关系,这些概率关系直接反映了变量间如何相互影响。

反映未知变量之间的条件依赖关系

在SLAM问题中,机器人的位置和地图的特征点都是未知变量,它们之间的关系通过传感器的测量数据(如激光雷达、摄像头等)被间接地捕获。因子图通过以下方式来反映这些未知变量之间的条件依赖关系:

  1. 变量节点:表示SLAM中的每一个未知元素,如每个时刻的机器人位置、地图上的每一个特征点。
  2. 因子节点:表示测量数据与这些未知变量之间的概率函数。例如,一个因子节点可以表示在给定机器人位置和方向的情况下,预期的激光雷达测量与实际测量之间的误差概率分布。

集中关注点在未知变量的推断或优化上

因子图将关注点放在了如何使用这些条件依赖关系来推断或优化未知变量上,而不是将测量数据本身作为解决问题的一部分。这是通过以下步骤实现的:

  1. 整合测量数据:测量数据作为参数整合到因子节点中,而不是作为独立的求解变量。这意味着我们不需要对测量数据本身进行优化或求解,而是使用这些数据作为已知的信息来帮助确定未知变量的最可能值。

  2. 优化或推断:使用图中的结构和因子定义,我们可以应用优化算法(如非线性最小二乘、梯度下降)或概率推断方法(如粒子滤波、贝叶斯推断)来估计这些未知变量的值。算法会根据因子图中定义的数学关系和输入的测量数据,寻找使得整个系统概率最大(或误差最小)的未知变量的配置。

通过这种方式,因子图能够有效地利用每个测量的信息,同时保持对未知变量间复杂关系的清晰表示,从而使得求解过程更为直接和高效。因此,因子图在SLAM及其他涉及大量变量和数据的领域中非常受欢迎。

在SLAM问题中,我们(通常)只考虑单个和成对的派系,导致一组未知Θ上的以下因子图表达式:

{{ $$ P(\Theta)\propto\prod_i\phi(\theta_i)\prod_{\{i,j\},i通常,势φ(θi)在未知θi∈Θ处编码先验或单一测量约束,而成对势ψij(θi,θj)与涉及两个未知θiand θj之间关系的测量或约束有关。注意,第二个产品是成对的团系{i,j},计算一次。方程(1)和(4)之间的等价性可以很容易地通过下列各式来建立 \[ \phi_0(x_0)\propto P(x_0)\\\psi_{(i-1)i}(x_{i-1},x_i)\propto P(x_i|x_{i-1},u_i)\\\psi_{i_kj_k}(x_{i_k},l_{j_k})\propto P(z_k|x_{i_k},l_{j_k}) \] 图1对应的完整SLAM问题的因子图表示。未知的姿态和地标分别对应于圆形和正方形的变量节点,而每个测量值对应于一个因子节点(填充的黑色圆圈)。

SLAM as a Markov Random Field

最后,用图形模型来表示SLAM问题的第三种方法是通过马尔可夫随机场,其中消除了因子节点本身。一个MRF的图是无向的,并且没有因子节点:它的邻接结构表示哪些变量是由一个共同的因子(度量或约束)连接起来的。在这个抽象级别上,形式(4)完全对应于一个成对的马尔可夫随机场[78]的表达式,因此是mrf和因子图是这里的等价表示。图1中的示例的MRF如图4所示。请注意,它看起来与图3非常相似,但是MRF是一个完全不同于因素图的结构(无向图和二部图)。

3.SAM as a Least Squares Problem

虽然前一节涉及建模,但我们现在讨论推理,即,给定我们可用的所有测量值,获得未知数集的最优估计。我们使用平滑而不是滤波,也就是说,我们想要恢复整个轨迹X={xi}和映射L={lj}的最大后验估计,给定测量值Z={zk}和控制输入U={ui}。让我们在Θ=(X,L)中收集X和L中的所有未知数。在上述假设下,我们通过最大化方程1中的联合概率P(X,L,Z)得到最大后验(MAP)估计。

\[ \begin{aligned}\Theta^*\triangleq\underset{\Theta}{\operatorname*{argmax}}P(X,L|Z)&=\quad\underset{\Theta}{\operatorname*{argmax}}P(X,L,Z)\\&=\quad\underset{\Theta}{\operatorname*{argmin}}-\log P(X,L,Z)\end{aligned} \]

这使我们通过(2)和(3)得出以下非线性最小二乘问题: \[ \Theta^*\triangleq\underset{\Theta}{\operatorname*{argmin}}\left\{\sum_{i=1}^M\|f_i(x_{i-1},u_i)-x_i\|_{\Lambda_i}^2+\sum_{k=1}^K\|h_k(x_{i_k},l_{j_k})-z_k\|_{\Sigma_k}^2\right\} \] 对于先前的P(x0),我们将假设x0是给定的,因此它被视为一个常数。这大大简化了本文档其余部分中的方程式。这是很常见的事情。在实践中:坐标系的原点是任意的,我们可以把x0固定在原点。该阐述很容易适应这种假设无效的情况。

在实践中,人们总是考虑公式(5)的线性化版本。如果过程模型和测量方程是非线性的,并且没有良好的线性化点,非线性优化方法,如高斯-牛顿迭代或莱文堡-马夸特算法将求解(5)的一系列线性近似,以接近最小[20]。这类似于[68,69,50]首创的针对SLAM的扩展卡尔曼滤波方法,但允许多次迭代以收敛,同时控制在哪个区域愿意信任线性假设(因此,这些方法通常被称为区域信任方法)。

我们现在将非线性最小二乘目标函数(5)中的所有项线性化。在最低点时,我们将假设一个好的线性化点是可用的,或者我们正在进行一个非线性优化方法的一次迭代。在任何一种情况下,我们都可以将(5)中的过程项线性化如下: \[ f_i(x_{i-1},u_i)-x_i\approx\left\{f_i(x_{i-1}^0,u_i)+F_i^{-1}\delta x_{i-1}\right\}-\left\{x_i^0+\delta x_i\right\}=\left\{F_i^{i-1}\delta x_{i-1}-\delta x_i\right\}-a_i\quad(6)\ \] \(\text{where }F_i^{i-1}\text{ is the Jacobian of }f_i(.)\text{ at the linearization point }x_{i-1}^0,\text{defined by}\) \[ F_i^{i-1}\triangleq\frac{\partial f_i(x_{i-1},u_i)}{\partial x_{i-1}}\Bigg|_{x_{i-1}^0}$​ \] \(a_i\stackrel{\Delta}{=}x_i^0-f_i(x_{i-1}^0,u_i)\)是测程预测误差(注意,\(u_i\)是这里是给出的,因此是常数)。(5)中的线性化测量项也同样得到, \[ h_k(x_{i_k},l_{j_k})-z_k\approx\left\{h_k(x_{i_k}^0,l_{j_k}^0)+H_k^{i_k}\delta x_{i_k}+J_k^{j_k}\delta l_{j_k}\right\}-z_k=\left\{H_k^{i_k}\delta x_{i_k}+J_k^{j_k}\delta l_{j_k}\right\}-c_k \]

\[ H_k^{i_k}\triangleq\frac{\partial h_k(x_{i_k},l_{j_k})}{\partial x_{i_k}}\bigg|_{(x_{i_k}^0,l_{j_k}^0)}J_k^{j_k}\triangleq\frac{\partial h_k(x_{i_k},l_{j_k})}{\partial l_{j_k}}\bigg|_{(x_{i_k}^0,l_{j_k}^0)}\ \]

\(c_k\triangleq z_k-h_k(x_{i_k}^0,l_{j_k}^0)\)是测量预测误差。

分别使用线性化的过程和测量模型(6)和(7),(5)开始 \[ \delta^*=\underset{\delta}{\operatorname*{argmin}}\left\{\sum_{i=1}^M\|F_i^{i-1}\delta x_{i-1}+G_i^i\delta x_i+a_i\|_{\Lambda_i}^2+\sum_{k=1}^K\|H_k^{i_k}\delta x_{i_k}+J_k^{j_k}\delta l_{j_k}-c_k\|_{\Sigma_k}^2\right\} \] 也就是说,我们得到了一个需要有效求解的δ中的线性最小二乘问题。为了避免处理\(\delta x_{i}\)的一种特殊方法,我们引入了矩阵\(G_i^i=-I_{d\times d},\),xi维数为d。

我们总是可以通过预乘以\(F_i^{i-1},G_i^i,\)\(a_i\)的每个项\(\Lambda_i^{-T/2},\)来消除\(\Lambda_i,\)(8),对于测量协方差矩阵Σk也是如此。对于标量测量,这仅仅意味着将每一项除以测量的标准差。下面我们假设已经这样做了,然后去掉马氏符号。

最后,将雅可比矩阵收集为矩阵a,并将向量纳入右侧(RHS)向量b,得到以下标准最小二乘问题, \[ \delta^*=\underset{\sim}{\operatorname*{argmin}}\left\|A\delta-b\right\|_2^2 \] 这是我们下面的起点。A可以变得非常大,但却非常稀疏,如图5所示。如果dx、dl和dz是状态、地标和测量值的尺寸,则A的大小为(Ndx Kdz)×(Ndx Mdl)。此外,A具有典型的块状结构,例如,具有M = 3、N = 2和K = 4: \[ \left.A=\left[\begin{array}{cccccc}G_1^1&&&&&\\F_2^1&G_2^2&&&&\\&F_3^2&G_3^3&&&\\H_1^1&&&J_1^1&&\\H_2^1&&&&J_2^2&\\&H_3^2&&J_3^1&&\\&&H_4^3&&J_4^2&\end{array}\right.\right],\quad b=\left[\begin{array}{c}a_1\\a_2\\a_3\\c_1\\c_2\\c_3\\c_4\end{array}\right] \] 上半部分描述了机器人的运动,下半部分是测量值。不同类型(和尺寸)的地标和/或测量值易于混合。

4. A Linear Algebra Perspective(线性代数透视)

cholesky 分解是什么

在本节中,我们将简要回顾 Cholesky 和 QR 因式分解及其在 (9) 中全秩线性最小二乘 (LS) 问题中的应用。(9) 中的全秩线性最小二乘 (LS) 问题的应用。这些材料已广为人知,给出它们主要是为了回顾线性代数算法与下一节的图论观点进行对比。论述紧跟[34],更深入的论述可参阅[34]。

原始问题是一个典型的最小二乘问题,目的是找到向量 \(\delta\),使得 \(A\delta - b\) 的欧几里得范数(即 \(L_2\) 范数)最小化:

\[ \delta^* = \underset{\delta}{\operatorname*{argmin}} \|A\delta - b\|_2^2 \]

最小化函数 \(\|A\delta - b\|_2^2\) 可以表示为:

\[ \|A\delta - b\|_2^2 = (A\delta - b)^T (A\delta - b) \]

展开此式:

\[ = (A\delta)^T (A\delta) - 2b^T (A\delta) + b^T b \]

\[ = \delta^T A^T A \delta - 2b^T A \delta + b^T b \]

为了找到这个函数的最小值,我们对 \(\delta\) 求导,并设其导数为零。由于这里是一个二次式,其导数可以用来找到极小值点:

\[ \frac{\partial}{\partial \delta} (\delta^T A^T A \delta - 2b^T A \delta + b^T b) = 0 \]

求导得:

\[ 2A^T A \delta - 2A^T b = 0 \]

简化此式,我们得到:

\[ A^T A \delta = A^T b \]

这就是所谓的正规方程,它提供了一个求解 \(\delta\) 的直接方法,即通过求解线性方程组 \(A^T A \delta = A^T b\) 来找到原始最小二乘问题的解。

这样,我们从最小化欧几里得范数的问题过渡到解线性方程的问题,两者实际上是等价的。解 \(A^T A \delta = A^T b\) 不仅找到了使得 \(A\delta\) 最接近 \(b\)\(\delta\),而且还以线性方程组的形式给出了直接的计算方法。这也是为什么在实际应用中(尤其是在数据拟合和机器学习中)我们经常需要计算 \(A^T A\)\(A^T b\)

Cholesky分解

对于全秩m×n矩阵a,具有m≥n,通过求解正规方程可以得到(9)的唯一LS解:\(A^TA\delta^*=A^Tb\)

这通常是通过信息矩阵I的Cholesky因子分解来完成的,定义和分解如下:\(\mathcal{I}\triangleq A^TA=R^TR\)

Cholesky三角分解 R是一个上三角n×n矩阵,使用Cholesky分解矩阵计算,即对称正定矩阵的LU分解的一个变体。对于密集矩阵,Cholesky分解需要n3/3的运算量。 \[ {\text{first }R^T}y=A^Tb\text{ and then }R\delta^*=y \] 在此之后,可以通过求解找到δ通过反向替代。

Cholesky分解需要n3/3的运算量。计算一半的对称ATA,需要mn2复杂度,整个算法,包括计算一半的对称ATA,需要(m+n/3)n2复杂度。

对于图1的例子,I和它的cholelesky三角形R在图5中与A一起显示。请注意,当A的列以典型的方式排序时,I的非常典型的块结构,例如,首先是轨迹X,然后映射L(下面我们将其称为XL排序) \[ \left.\mathcal{I}=\left[\begin{array}{cc}A_X^TA_X&\mathcal{I}_{XL}\\\mathcal{I}_{XL}^T&A_L^TA_L\end{array}\right.\right] \] \(\mathcal{I}_{XL}\triangleq A_X^TA_L\)对机器人状态X和映射L之间的相关性进行编码,对角线块为带对角线。

避免计算平方根的Cholesky分解的一个变体是LDL分解,它计算一个下三角矩阵L和一个对角矩阵D \[ \mathcal{I}=R^TR=LDL^T \] 另一种更精确和数值稳定的Cholesky分解的方法是通过QR分解而不计算信息矩阵I。相反,我们计算A本身的QR分解及其相应的RHS: \[ \left.Q^TA=\left[\begin{array}{c}R\\0\end{array}\right.\right]\quad Q^Tb=\left[\begin{array}{c}d\\e\end{array}\right] \] 这里Q是一个m×m正交矩阵,R是上三角的Cholesky三角。分解密集矩阵A的首选方法是逐列计算R,从左到右进行。对于每一列j,对角线下的所有非零元素都通过将左边的A与Householder反射矩阵Hj相乘而归零。经过n次迭代后,A被完全分解。 \[ \left.H_n..H_2H_1A=Q^TA=\left|\begin{array}{c}R\\0\end{array}\right.\right| \] 正交矩阵Q通常不形成:相反,转换后的RHS \(Q^Tb\)是通过将b作为A的附加一列来计算的。因为Q因子是正交的,我们有: \[ \left\|A\delta-b\right\|_2^2=\left\|Q^TA\delta-Q^Tb\right\|_2^2=\left\|R\delta-d\right\|_2^2+\left\|e\right\|_2^2 \] 显然,||e||将是最小二乘残差,通过求解平方系统可以得到LS解δ∗:Rδ = d

通过反向替代。QR的成本主要由Householder反射的成本决定,即2(m−n/3)n2

比较QR和Cholesky分解,我们发现两种算法都需要O(mn2)运算(m>n),但QR分解要慢2倍。虽然这些数字只对密集矩阵有效,但我们已经看到,在实践中,LDL和Cholesky因子分解在稀疏问题上也远远优于QR因子分解,而不仅仅是一个常数因子。

5.A Graphical Model Perspective(图像模型透视图)

5.1 Matrices ⇔ Graphs (矩阵 图)

从上面的说明中,现在可以很容易地认识到,测量雅可比矩阵A是与SLAM相关的因子图的矩阵。我们可以从两个层面上理解这个说法。首先,A的每个块对应于最小二乘准则(8)中的一项,要么是一个地标测量,要么是一个测程约束,并且每一个块行对应于因子图中的一个因子。在每个块行中,稀疏性模式表示哪些未知的姿态和/或地标与该因素相连。因此,A的块结构完全对应于与SAM相关的因子图的邻接矩阵。

其次,在标量水平上,每一行Ai在A(见图5)对应于稀疏矩阵最小二乘准则(9)中的一个标量项\(\|A_i\delta-b_i\|_2^2\),如 \[ \left\|A\delta-b\right\|_2^2=\sum_i\left\|A_i\delta-b_i\right\|_2^2 \] 因此,这定义了一个精细结构的因图,通过 \[ P(\delta)\propto\exp-\frac12\left\|A\delta-b\right\|_2^2=\prod_i\exp-\frac12\left\|A_i\delta-b_i\right\|_2^2 \] 重要的是要认识到,在这个更精细的观点中,SLAM问题的块结构被丢弃了,并且它是这个图被通用的线性代数方法所检验。通过使用块结构来代替,我们将能够做得更好。

如之前在[75]以及其他文献中提到的,信息矩阵Ⅰ=ATA是与SLAM问题关联的马尔可夫随机场(MRF)的矩阵表达形式。同时,从块级别来看,ATA的稀疏模式完全等同于相应MRF的邻接矩阵。方程5中的目标函通过哈默斯利-克利福德定理[77]对应于成对马尔可夫随机场(MRF)[77,78]。而在这个MRF中,节点对应于机器人的状态和地标。链接代表里程计或地标测量。

在文献[63, 75]中,采用马尔可夫随机场(MRF)图的视角来揭示SLAM滤波版本中固有的相关性结构。研究表明,当将过去的轨迹 \(𝑋_{1:𝑀−1}\)边缘化时,信息矩阵不可避免地变得完全密集。因此,这些方法的重点是选择性地移除链接以减少滤波器的计算成本,并取得了显著的成功。相比之下,本文考虑的是与平滑信息矩阵\(I\)关联的MRF,该信息矩阵不会变得密集,因为过去的状态从未被边缘化。

5.2 因式分解⇔变量消除

剩下的一个问题是平方根信息矩阵R对应于什么图?请记住,R是在第4节中对I或A进行因素分解的结果。Cholesky或QR因子分解最常被用作“黑盒”算法,但实际上它们与最近在图形模型[11]中开发的推理方法相似。从下面可以看到,R与结树基本对应,从图形模型的推理中知道,最近也在SLAM中应用

两种因式分解方法,QR和Cholesky(或LDL),都基于变量消除算法[4,11]。这两种方法的区别在于QR通过消除因子图中的变量节点得到\[A=QR\],而Cholesky或LDL则从马尔可夫随机场(MRF)开始,因此得到\[A^TA = R^TR\]。这两种方法都是从\[\delta_1\]开始,一次消除一个变量,\[\delta_1\]对应于\[A\]\[\mathcal{I}\]的最左列。消除的结果是\[\delta_1\]现在表示为所有其他未知数\[\delta_{j>1}\]的线性组合,系数位于对应的\[R_1\]行中。然而,在此过程中,新的依赖性被引入到所有与\[\delta_1\]相连的变量之间,这导致图中添加了边。然后以类似方式处理下一个变量,直到所有变量都被消除。这正是图形模型推理中熟悉的公理化和三角化过程。消除所有变量的结果是一个有向的、三角化(弦图)的图,在我们的示例中显示在图6中。

这段话描述的是在处理线性系统或矩阵分解时,如何使用变量消除算法,具体包括QR分解和Cholesky分解(或LDL分解)。以下是每部分的具体解释:

  1. 变量消除算法: 这是一种算法框架,用于通过逐步简化变量来处理数学问题。在矩阵分解的上下文中,这意味着选择某个变量并重新组织方程,使得这个变量可以通过其它变量线性表示。

  2. QR分解和Cholesky分解的对比:

    • QR分解(\[A = QR\]):这是通过在因子图中消除变量节点实现的。在此过程中,矩阵\[A\]被分解为一个正交矩阵\[Q\]和一个上三角矩阵\[R\]。QR分解特别适用于求解具有线性独立列的线性系统。
    • Cholesky分解(\[A^TA = R^TR\]):这种方法以马尔可夫随机场(MRF)为起点。在数学上,这涉及到将矩阵\[A\]的转置与\[A\]本身相乘,得到的结果是一个对称正定矩阵,然后对这个矩阵进行因式分解以得到一个上三角矩阵\[R\]。此方法适用于正定矩阵的分解。
  3. 变量消除过程:

    • \[\delta_1\]开始,它是矩阵\[A\]\[\mathcal{I}\]的最左列所对应的变量。这个变量被消除或解算出来,并用其它变量(\[\delta_{j>1}\])的线性组合表示。
    • 在消除\[\delta_1\]的过程中,会在图中介入新的依赖性,这些依赖性在与\[\delta_1\]相连的所有变量之间形成新的边。这是因为消除一个变量通常会使其相关的其他变量之间产生直接的联系,这些联系在原始图中可能并不存在。
  4. 公理化和三角化过程:

    • 这是图形模型推理中的常见步骤。消除变量的过程不仅涉及代数运算,还涉及图的结构变化,比如添加新的边来填充环路,使得图变成有向且无环的。
  5. 结果:

    • 完成所有变量的消除后,结果是一个有向的、三角化(或称为弦图)的图。在此上下文中,这意味着所有的变量都通过某种方式关联,并且关联结构支持高效的图形模型推理。

这些概念通常用于数值分析、统计建模和计算机视觉等领域,用于优化算法的效率并减少计算复杂性。

一旦获得了弦图 \((R!)\),我们就可以得到 \(R\) 的消元树,它被定义为消元后弦图的深度优先生成树,并且在反向替代阶段中用于说明计算流程。为了说明这一点,图 7 展示了使用先消除地标再消除姿态的著名启发式方法所获得的弦图(我们将其称为 \(LX\) 排序)。相应的消元树在图 8 中展示。树的根对应于最后一个被消除的变量 \(\delta_n\),它是在反向替代(方程 13)中首先被计算的。然后,计算沿树向下进行,尽管这通常按逆列顺序完成,但在不相连的子树中的变量可以按任何顺序计算。实际上,如果只对某些变量感兴趣,那么没有必要计算不包含这些变量的任何子树。

然而,分析并没有就此结束。\(R\) 的图具有一个可以完全封装在称为团树(clique tree)的根树数据结构中的团结构 [64, 4],在人工智能文献中也称为连接树(junction tree)[11]。例如,图 1 问题上的 \(LX\) 排序的团树在图 9 中展示。对应关系几乎是一对一的:每个 \(R\) 准确对应一个团树,反之亦然,仅需在团内进行列重排。团树也是多前端 QR 方法 [56] 的基础,我们在下面的模拟中也对此进行了评估。在多前端 QR 因式分解中,计算从树的叶子到根进行,以分解一个稀疏矩阵,然后从根到叶子进行反向替代步骤。关于平方根信息矩阵和团树之间关系的完整处理超出了当前论文的范围,但在其他工作中,我们已经在分布式推断的新算法中使用了团树结构 [19]。

5.3提高性能⇔减少填料

良好性能的最重要的因素是变量的顺序。不同的变量顺序可以显著地产生或多或少的填充,定义为在因子分解过程中添加到图中的边的量。由于每条边对应于Choleseky三角形R中的非零,计算R和反向替换的成本在很大程度上取决于发生的填充量。不幸的是,找到一个最优排序是np完备的。 发现近似最优排序的算法是稀疏线性代数研究的一个活跃领域。一种常用的处理中型问题的方法是联合[1],另一种基于图论的常用方法是广义嵌套解剖[53,52]。

假设最优排序通常是遥不可及的,启发式或领域知识可以比通用算法做得更好。一个简单的想法是使用一个标准的方法如colamd,但它工作块的稀疏模式而不是通过原来的测量雅可比a正如我们上面提到的,块结构是真正的知识SLAM问题,不能访问colamd或任何其他近似排序算法。 虽然对colamd性能的影响可以忽略不计,但我们发现,让它在SLAM MRF上工作,而不是直接在稀疏矩阵上工作,有时可以提高2到100倍,15是一个很好的经验法则。

请注意,在某些情况下,任何排序都会导致相同的大填充。最坏的情况是一个完全连接的二部2MRF:每个地标都可以从每个位置看到。在这种情况下,消除任何变量将完全连接所有变量的另一边,之后集团树的结构是完全已知的:如果首先选择一个姿势,根将整个地图,和所有姿态将计算一旦地图。反之亦然,如果选择了一个地标,那么轨迹将是团树根团,计算将通过(昂贵的)轨迹优化进行,然后是(非常便宜的)地标计算。最重要的是,这两种情况构成了“舒尔补体”技巧的基础,它在结构中从运动应用[76,40]中闻名,也在GraphSLAM [74]中使用。

然而,上面概述的最坏情况是机器人技术中的一个例外:传感器的范围有限,并且被墙壁、物体、建筑物等遮挡。这在大规模映射应用程序中尤其如此,它本质上意味着MRF通常将是稀疏连接的,即使它是一个大的连接组件。

这一段文本详细解释了如何利用QR和Cholesky(或LDL)因子分解方法来处理矩阵,以及这些方法与图形模型中的推理如何相互关联。下面我将逐步解释这些内容:

  1. 因子分解方法:
    • QR 和 Cholesky(或LDL)因子分解:这两种方法都基于变量消除算法。QR方法通过从因子图中消除变量节点来获得 \[A = QR\],而Cholesky或LDL方法从马尔可夫随机场(MRF)出发,得到 \[A^T A = R^T R\]
    • 变量消除过程:从 \[\delta_1\] 开始,逐一消除变量,\[\delta_1\] 是矩阵 \[A\]\[\mathcal{I}\] 的最左列。消除的结果是,\[\delta_1\] 被表示为所有其他未知数 \[\delta_{j>1}\] 的线性组合,系数存放在 \[R_1\] 行中。这个过程引入了新的依赖性,导致图中增加了边。
  2. 图的形成与消元树:
    • 形成弦图:通过上述变量消除,最终得到一个有向的、三角化(弦图)的图。
    • 消元树:一旦得到弦图,可以定义一个消元树,这是消元后弦图的深度优先生成树。这个树在反向替代阶段中非常有用,用于指导计算流程。
  3. 团树结构:
    • 团树(clique tree)或连接树(junction tree)\[R\] 的图具有可以完全封装在团树数据结构中的团结构。团树几乎与 \[R\] 一一对应,每个 \[R\] 精确对应一个团树,反之亦然,只需在团内进行列重排。
    • 多前端QR方法:在多前端QR因式分解中,计算从树的叶子到根进行,以分解一个稀疏矩阵,然后从根到叶子进行反向替代步骤。这种方法基于团树结构进行计算。
  4. 应用与扩展:
    • 在其他工作中,已经在分布式推断的新算法中使用了团树结构。这显示了此类数学工具不仅限于单一应用,而是可以广泛应用于各种复杂问题的求解中。

这段描述不仅阐述了两种主要的因子分解方法的技术细节,还链接了这些数学技术和图形模型理论的关系,尤其是在如何使用图结构来优化计算和推理过程中。

6 Square Root SAM

在本节中,我们将获取我们从上面所知道的所有内容,并说明三个简单的√SAM变体,这取决于它们是以批处理还是增量模式运行,以及是否涉及非线性。

6.1 Batch√SAM

一个批量版本的平方根信息平滑和映射是很简单的,也是解决一个大的、稀疏的最小二乘问题的一种通用的标准方法:

\[ \begin{aligned}&\frac{\textbf{Algorithm 1 Batch }\sqrt{\text{SAM}}}{\text{I. Build the measurement Jacobian }A\mathrm{~and~the~RHS~}b\text{ as explained in Section 3.}}\\&2.\text{ Find a good column odering }p,\text{ and reorder }A_p\xleftarrow{p}A\\&3.\text{ Solve }\delta_p^*=\text{ argmin }_\delta\left\|A_p\delta_p-b\right\|_2^2\text{ using ether the Cholesky or QR factorization method}\\&\text{from Section 4}\\&4.\text{ Recover the optimal solution by }\delta\xleftarrow{r}\delta_p,\text{with }r=p^{-1}\end{aligned} \]

  1. 建立第3节所述的雅可比值A和RHS b
  2. 找到一个好的列排序p,并重新排序App←a
  3. 使用第4节中的Cholesky或QR分解方法 解决 δ∗p=argminδkApδp−bk22
  4. 用δr←δp恢复最优解,用r=p−1

在测试中,我们获得了稀疏LDL分解[13]的最佳性能,如上所述,它是计算I = LDLT的一个变体,D为对角矩阵,L是对角上的三角矩阵。 同样的算法也适用于非线性情况,但被非线性优化器简单地调用。但是,请注意,排序只需要计算一次。

其他

消元树

消元树是一种用来表示消除过程中变量依赖关系的树状结构,它帮助指导了矩阵分解过程中的计算流程,特别是在反向替代阶段。这里用一个具体的例子来说明这一点:

假设我们有一个系统的矩阵 \[A\],我们想要进行Cholesky分解(假设\[A\]是对称正定的)。我们首先选择一个变量来消除,这通常是基于某种启发式方法,比如最小填充法或最小度法。

考虑以下矩阵:

\[ A = \begin{bmatrix} 4 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix} \]

第一步:消除第一个变量

我们首先消除第一个变量,这会影响与它直接相连的变量(这里是第二个变量)。消除第一个变量后,更新矩阵(做适当的行和列操作)可能会导致非零元素填充原本为零的位置,即产生“填充”。

弦图的形成

在这个过程中,如果我们将每个变量和与它相邻的变量视为图的一个节点和边,消除一个变量相当于在图中移除一个节点及其相关的边,并可能在未被消除的相邻节点间添加新的边(如果它们之前不直接相连)。这个操作最终形成了一个弦图,即每个循环都至少有一个弦(非环的边)的图。

构建消元树

在构建弦图的基础上,消元树被定义为消除过程中的依赖关系图。例如:

  1. 将第一个被消除的变量作为树的根节点。
  2. 每次消除一个变量后,将受到影响的变量作为当前消除变量的子节点添加到树中。

对于上面的矩阵,消元树可能如下:

  • 根节点为第一个被消除的变量。
  • 第二个被影响的变量成为第一个节点的子节点。
  • 依次类推,直到所有变量都被放置在树中。

在反向替代阶段的作用

在进行例如Cholesky分解后,求解线性系统 \[Ax = b\] 需要先进行前代(使用 \[L\] 矩阵)和后代(使用 \[L^T\] 矩阵)。在后代过程中,消元树指导了计算的顺序:从树的叶子(最先消除的变量)开始,逐步向根节点(最后消除的变量)进行,确保在计算当前变量值之前,所有依赖于它的变量值已经被计算。

通过这种方式,消元树不仅帮助优化计算过程,还确保了在计算过程中遵循正确的依赖关系,有效地利用了矩阵的稀疏性质。

滤波器与平滑器

在SLAM(Simultaneous Localization and Mapping,即同时定位与建图)问题中,信息矩阵 (I) 是一个关键的概念,用于表示机器人的状态(如位置和方向)以及环境中地标的不确定性。信息矩阵是协方差矩阵的逆,其中每个元素反映了状态变量之间的条件依赖性。

滤波器与平滑器的区别

在SLAM中,主要有两种处理信息的方法:滤波和平滑。

  1. 滤波
    • 在滤波过程中,为了保持计算的实时性和可管理性,通常只会保留当前状态的信息,而把旧的状态信息“边缘化”(即从当前考虑的状态集中移除)。这种边缘化过程会引入额外的非零元素到信息矩阵中,因为边缘化旧状态意味着需要将这些状态的信息“分摊”到剩余的状态上。结果是信息矩阵变得密集,即使原来很多状态之间是条件独立的。
    • 为了处理这种密集化,研究者通常会采用各种技术(如稀疏化技术)来移除信息矩阵中的某些链接,以减少计算负担。
  2. 平滑
    • 平滑处理,相比于滤波,不会边缘化任何过去的状态。相反,它保留了所有历史和当前的状态信息,因此可以在整个时间线上提供更精确的状态估计。
    • 在平滑方法中,由于不进行边缘化操作,所以不会引入额外的非零元素导致信息矩阵密集化。这意味着信息矩阵能够保持其原始的稀疏结构,其中非零元素主要集中在时间线上相邻状态之间的依赖关系。这种稀疏性是因为在大多数SLAM场景中,只有时间上相邻的状态或空间上邻近的地标才具有直接的相互依赖。

结论

因此,在与平滑信息矩阵 (I) 关联的MRF中,信息矩阵保持稀疏,主要是因为平滑处理不涉及边缘化过去的状态。这样,每个状态只与其直接相关的少数几个其他状态或地标相关联,而不是与整个历史状态集相关联,从而避免了矩阵的密集化。这种方法不仅保持了计算效率,而且通常能提供比滤波更为精确的状态估计。

滤波是一种在线处理方式,即它只处理当前和过去的信息,通常用于实时系统中。在滤波过程中,随着新数据的持续加入,为了维持计算的可管理性,旧的数据或状态往往需要被边缘化(移除),这导致信息矩阵趋向于变得密集。

平滑则是一种离线处理方式,它考虑全部的数据(包括过去和未来的数据)。在平滑处理中,由于不需要边缘化任何状态(因为所有状态对计算都是可用的),所以关联的信息矩阵维持稀疏性,这是因为不存在需要通过边缘化引入新依赖的需求。

滤波

当然可以,让我们通过一个具体的例子来解释为什么在SLAM中进行滤波操作时,边缘化过去的状态会导致信息矩阵变得密集。

示例情景

假设一个机器人在一个简单的环境中沿直线移动,每移动一段距离记录一次位置,形成状态序列 (X_1, X_2, X_3, , X_t)。每个状态都与其前后状态有直接的连接(例如通过里程计测量),并可能与环境中的某些地标有关联。

在开始,信息矩阵 (I) 是稀疏的,因为只有时间上相邻的状态或与某个地标直接相关的状态之间存在连接。例如,如果(X_1)和(X_2)通过里程计相关联,那么在信息矩阵中,(I_{1,1}),(I_{1,2}),和(I_{2,2})是非零的。

边缘化过程

当机器人继续前进并且需要在滤波器中处理新的状态时(比如(X_{t+1})),为了保持计算的可管理性,较早的状态(比如(X_1)) 可能会被边缘化掉。边缘化(X_1)意味着我们需要从当前激活的状态集中移除(X_1)的影响。

  • 原本(X_1)直接影响(X_2),在边缘化(X_1)后,我们必须将(X_1)对(X_2)的影响“转移”给与(X_2)相关的其他状态,比如(X_3)。这意味着原先(X_1)与(X_3)可能没有直接关系,但在边缘化(X_1)后,(X_2)与(X_3)之间的关系需要加强,以反映从(X_1)到(X_3)的间接影响。

结果

这个转移过程在信息矩阵中表现为原先为零的元素变成非零,因为(X_2)和(X_3)现在有了新的依赖关系。随着越来越多状态的边缘化,这种间接依赖会累积,导致信息矩阵逐渐从稀疏变为密集。

几何直观

可以想象,每个状态不仅与其直接的邻居有联系,还因为边缘化而与更多其他状态产生了联系。随着时间的推移,几乎每个状态都以某种方式间接地与其他所有状态连接,从而增加了信息矩阵中的非零元素数量,使矩阵变得密集。

这种密集化不仅增加了计算的复杂度,还可能影响滤波器的性能,因此在实际应用中,开发者会寻找各种方法(例如使用稀疏近似或选择性地保留关键状态)来尽可能减少这种影响。

平滑

我们可以用一个具体的数学例子来解释SLAM问题中平滑方法的信息矩阵是如何维持其稀疏性的。在这个例子中,我们考虑一个小规模的机器人轨迹,其中机器人在四个时间点上测量其位置。我们将看到在不进行边缘化的情况下,信息矩阵如何保持其结构。

假设有四个状态 (X_1, X_2, X_3, X_4),它们代表机器人在四个连续时间点的位置。在SLAM中,通常只有连续的状态之间存在直接的测量依赖关系。

在平滑方法中,信息矩阵 (I) 初始化时将反映状态之间的这些直接依赖关系。例如,如果我们假设每个状态只与其前后状态有关系,并且每个这样的依赖对信息矩阵贡献 (1) 的信息值,那么信息矩阵可能看起来像这样:

\[ I = \begin{bmatrix} 1 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 1 \end{bmatrix} \] 这里,对角线上的元素 (2) 表示一个状态依赖于其前后状态的信息总和(除了边界条件,即 (X_1) 和 (X_4))。非对角线的 (-1) 表示状态 (X_i) 与 (X_{i+1}) 和 (X_{i-1}) 之间的信息依赖。

在平滑方法中,由于不进行边缘化操作,我们保留所有历史状态的数据。这意味着在任何给定的时刻,我们都可以访问到所有历史和未来的状态信息,而不是只有当前和未来的状态。因此,信息矩阵的非零结构不会改变——它始终保持与初始时刻相同的稀疏模式。

这种方法的优点是,信息矩阵维持了其稀疏性,这大大减少了计算的复杂性和内存需求。例如,当使用稀疏矩阵技术处理大型SLAM问题时,这一特性非常重要。每个状态只与直接相邻的状态有关系,没有因边缘化而引入的额外复杂依赖关系。

这样的稀疏结构是高效处理大规模SLAM问题的关键,因为它减少了必须处理的数据量,并允许使用优化的数值方法来处理信息矩阵。

为什么在平滑中仍然会有变量消除

尽管在平滑处理中,信息矩阵维持其稀疏性质,我们仍可能需要执行某种形式的变量消除,这主要是为了效率和算法实现上的考虑。例如,在求解线性系统时(如 \[Ax = b\]),通常会采用一些数学技术来简化问题,比如通过Cholesky分解等。在这个过程中,创建弦图和消元树是为了优化这些数学操作:

  1. 弦图的创建:在进行Cholesky分解等操作前,弦图能帮助我们识别并利用信息矩阵的稀疏结构,以减少计算过程中的填充(非零元素的添加)。

  2. 消元树的建立:即便信息矩阵是稀疏的,构建消元树可以进一步提高解线性系统的效率。消元树能够指导计算顺序,确保依赖关系得到正确处理,并优化内存访问模式。

因此,尽管平滑信息矩阵在理论上不需要通过边缘化来维护其稀疏性,但在实际计算和算法实现中,变量消除和相关的图形表示(如弦图和消元树)仍然是优化性能和保证算法效率的重要工具。

QR分解

QR分解是将矩阵 ( A ) 分解成一个正交矩阵 ( Q ) 和一个上三角矩阵 ( R ) 的方法。如果 ( A ) 是一个 ( m n ) 矩阵且 ( m n ),那么分解形式通常写为:

\[ A = QR \]

其中,( Q ) 是 ( m m ) 正交矩阵(( Q^T Q = QQ^T = I )),( R ) 是 ( m n ) 的上三角矩阵,其中除了前 ( n ) 行可能有非零元素外,其余部分全为零。在许多实际情况中,我们更关心 ( R ) 的前 ( n n ) 部分,因为它是一个完全的上三角矩阵。

为什么 ( R ) 下面有一个0

在 ( Q^T A = ) 的表达式中,( R ) 是 ( n n ) 的上三角矩阵,而下面的0表示在 ( Q^T ) 作用后,剩下 ( m-n ) 行都被消减到了零。这是因为 ( Q ) 通过正交变换将 ( A ) 的所有向量投影到一个由 ( A ) 的列空间张成的正交子空间。

如何进行QR分解

QR分解可以通过多种方式进行,包括Householder反射、Gram-Schmidt正交化或Givens旋转等。其中,Householder反射是最常用的方法,因为它具有很好的数值稳定性。

Householder反射:

  1. 选取一个列向量,构造一个反射矩阵 ( H ),使得除了第一项外,该列的其余项变为零。
  2. 重复应用上述过程,每次将 ( A ) 的一个列向量(逐步从左到右)变为期望的上三角形式。
  3. 通过连续左乘不同的Householder矩阵(( H_1, H_2, , H_n )),我们最终得到 ( Q^T A = )。

QR分解在最小二乘中的应用

在最小二乘问题中,我们通常要最小化 ( |A- b|_2^2 )。通过QR分解,我们可以转化这个问题为:

\[ \|Q^T(A\delta - b)\|_2^2 = \|R\delta - Q^Tb\|_2^2 \]

由于 ( Q ) 是正交的,所以转换不会改变二范数。因此,问题简化为最小化 ( |R- d|_2^2 )(其中 ( d ) 是 ( Q^Tb ) 的前 ( n ) 个元素)加上 ( |e|_2^2 )(其中 ( e ) 是 ( Q^Tb ) 的剩余 ( m-n ) 个元素,代表方程的残差)。

由于 ( R ) 是上三角的,这个问题可以通过反向替换直接求解 ( R= d ),从而找到最小二乘解 ( ^* )。

在说明为什么在 ( A = QR ) 分解中,矩阵 ( R ) 除了前 ( n ) 行可能有非零元素外,其余部分全部为 0 之前,有必要澄清上面的描述有些不够精确。对于标准的 ( A = QR ) 分解,特别是当 ( A ) 是 ( m n ) 矩阵且 ( m n ) 时,我们应该更明确地表示这个关系。

正确的 ( R ) 矩阵结构

在 ( A = QR ) 分解中,( R ) 实际上是一个 ( n n ) 的上三角矩阵,而不是 ( m n ) 的矩阵。这是因为 ( Q ) 是一个 ( m m ) 的正交矩阵,而 ( R ) 是一个 ( m n ) 的矩阵,其中只有顶部的 ( n n ) 部分是上三角形的,其余 ( m-n ) 行则全部为零。这样,( QR ) 的乘积仍然能重构出原始的 ( m n ) 矩阵 ( A )。

为什么 ( R ) 的下部为零

  1. 正交变换的作用:在 ( QR ) 分解中,( Q ) 的作用是通过正交变换(如Householder变换或Givens旋转)逐步将 ( A ) 的列向量变为上三角形式。每次变换主要集中在将当前处理的列下方的元素变为零。

  2. 维度限制:当处理到 ( A ) 的第 ( n ) 列时,如果 ( m > n ),则在 ( n ) 行之下没有更多的列元素可以继续处理。因此,在完成所有列的变换后,( R ) 矩阵中 ( n ) 行以下的部分将会保持为零。这是因为这部分 ( Q ) 中的变换已经将 ( A ) 的所有相关信息“推”到了其前 ( n ) 行中的上三角形状。

因此,( R ) 矩阵的结构确保了 ( QR ) 分解的过程只影响 ( A ) 的前 ( n ) 行,而其余 ( m-n ) 行为零,确保 ( QR ) 乘积可以精确重构原始矩阵 ( A )。这种结构特点是 ( QR ) 分解的一个重要特性,使得它在解决线性最小二乘问题和其他矩阵相关计算中非常有效。

信息矩阵I

在机器人定位和地图构建的问题中,通常涉及到估计机器人在环境中的位置(轨迹 (X))以及环境本身的结构(地图 (L))。在这种情况下,信息矩阵 (),也被称为费舍尔信息矩阵,是用来表示这些未知量的不确定性和相互依赖关系的关键工具。

信息矩阵 () 是由雅可比矩阵 (A)(或系统的设计矩阵)通过 (A^T A) 计算得到的。这个矩阵可以从以下几个方面描述机器人轨迹和地图内部的信息:

  1. 量化不确定性:
    • 对角线元素:信息矩阵的对角线元素量化了每个参数(位置或地图特征)的信息量或确信度。对角线元素的值越大,表示对应参数的不确定性越小,我们对该参数的估计越精确。
    • 非对角线元素:非对角线元素描述了参数之间的依赖关系或相关性。例如,如果位置 (x_i) 和位置 (x_j)(或地图特征)在矩阵中的对应元素非零,这表明这两个位置(或特征)的估计值相互依赖。
  2. 编码相关性:
    • 分块 ({XL}) 和 ({XL}^T) 特别重要,因为它们编码了轨迹 (X) 和地图 (L) 之间的相关性。这意味着机器人在特定位置的估计如何依赖于地图的特定特征,反之亦然。这种相关性对于同时定位和地图构建(SLAM)至关重要,因为它帮助系统同时解决导航和地图创建的问题。
  3. 优化和解算:
    • 在实际应用中,基于信息矩阵的结构,可以采用更有效的数值方法来解决估计问题,如使用稀疏矩阵技术处理大规模SLAM问题。信息矩阵的稀疏性表明很多参数之间的依赖性较低,这可以大大减少计算量。

总之,信息矩阵 () 不仅提供了对机器人轨迹和地图每一部分的内部信息的详尽描述,而且还通过其结构揭示了不同部分之间的相互作用和依赖性,从而在理论和实践中都是解决SLAM问题的核心部分。

这段话描述的是矩阵 (),它代表了信息矩阵,在这个具体例子中,它是用于机器人状态估计的问题。矩阵 () 是通过矩阵 (A) 的转置乘以 (A) 计算得到的,这里 (A) 被分为两部分:对应于机器人轨迹 (X) 的 (A_X) 和对应于地图 (L) 的 (A_L)。

矩阵 () 的结构表明它分为四个部分:

  1. (A_X^T A_X):这是左上角的块,代表轨迹 (X) 内部的信息。
  2. (A_L^T A_L):这是右下角的块,代表地图 (L) 内部的信息。
  3. ({XL}) 和 ({XL}^T):这两个块分别位于非对角线的位置,(_{XL} A_X^T A_L) 和其转置表示轨迹 (X) 和地图 (L) 之间的相关性。

\[ \mathcal{I}_{XL} \triangleq A_X^T A_L \]

表示的是 (A_X) 和 (A_L) 之间的交叉信息,它编码了机器人状态 (X) 和地图 (L) 之间的相关性。对角线块 (A_X^T A_X) 和 (A_L^T A_L) 表明它们是“带对角线”的,这通常意味着这些块主要关注各自的内部结构和相互作用,而不是外部的。

关于AX-b 中A和b的物理意义

在SLAM(Simultaneous Localization and Mapping)问题中,机器人的动态与环境交互通过一系列的状态转移和观测更新来描述。矩阵 ( A ) 在这个问题中的作用,是将这些动态和交互的关系转化为数学形式,特别是以线性化的方式来近似非线性系统的行为。下面详细解释这一点:

机器人的动态

机器人的动态指的是机器人如何根据其控制指令(如速度和方向)从一个状态(位置和方向)移动到另一个状态。这些动态通常可以用一个非线性函数 ( f ) 表示:

\[ x_i = f(x_{i-1}, u_i) \]

其中,( x_i ) 和 ( x_{i-1} ) 分别是连续两个时间步的状态,( u_i ) 是在时间 ( i ) 应用的控制指令。在SLAM问题的线性化处理中,我们对函数 ( f ) 在某个点进行泰勒展开,通常是在预测状态附近,得到:

\[ x_i \approx F_i^{i-1} x_{i-1} + G_i^i u_i \]

其中,( F_i^{i-1} ) 和 ( G_i^i ) 是偏导数矩阵,代表状态和控制输入对新状态的线性影响。这一部分在矩阵 ( A ) 中以块的形式出现,描绘了状态之间的转移关系。

环境的观测

环境的观测描述了机器人如何通过其传感器(如摄像头、激光雷达)观察到周围环境,特别是地标的位置。这些观测同样可以表示为一个非线性关系 ( h ):

\[ z_k = h(x_{i_k}, l_{j_k}) \]

其中,( z_k ) 是第 ( k ) 次观测,( x_{i_k} ) 是观测时的机器人状态,( l_{j_k} ) 是被观测的地标。同样地,通过在预测观测点附近进行泰勒展开,我们得到观测的线性化模型:

\[ z_k \approx H_k^{i_k} x_{i_k} + J_k^{j_k} l_{j_k} \]

这里的 ( H_k^{i_k} ) 和 ( J_k^{j_k} ) 代表观测模型的线性化矩阵,反映了状态和地标位置对观测结果的线性影响。这部分在矩阵 ( A ) 中也以块的形式呈现。

总结

因此,矩阵 ( A ) 通过其块状结构和各个元素,捕捉了机器人如何移动(动态)和如何观察周围环境(观测)的数学描述。这不仅包括了机器人状态的转移关系,还包括了状态与环境地标之间的相互作用。通过优化问题最小化 ( |A- b|_2^2 ),我们实际上是在调整对机器人状态和地图的估计,使得这些线性化的预测尽可能贴近实际观测,从而实现同时定位与地图构建的目的。 在SLAM(Simultaneous Localization and Mapping)问题中,矩阵 ( A ) 和向量 ( b ) 代表了机器人定位与地图构建问题的线性化模型。这里详细解释矩阵 ( A ) 的意义和重要性:

矩阵 ( A ) 的意义:

  1. 物理意义:矩阵 ( A ) 表示系统状态变量之间的关系,以及状态变量与观测之间的线性关系。具体来说:
    • 运动模型部分:矩阵 ( F_i^{i-1} ) 和 ( G_i^i ) 表示从状态 ( x_{i-1} ) 到状态 ( x_i ) 的线性化运动模型。这些块表示了机器人如何从一个位置(和/或方向)移动到下一个位置,通常基于机器人的控制指令如转向和前进。
    • 观测模型部分:矩阵 ( H_k^{i_k} ) 和 ( J_k^{j_k} ) 表示状态 ( x_{i_k} ) 和地标 ( l_{j_k} ) 如何影响第 ( k ) 次观测。这反映了如何从当前状态和地标位置预测观测值(例如,从摄像头或激光雷达)。
  2. 计算意义:在最优化问题中,我们通常求解最小化 (|A- b|_2^2),这是一个最小二乘问题,意味着我们在寻找最佳的状态变量更新 ( ),使得线性化模型的预测尽可能接近实际观测和运动指令的结果。

为什么 ( A ) 重要:

  • ( A ) 矩阵捕捉了机器人动态和环境的整体结构,这是理解和解决SLAM问题的核心。正确构建和更新这个矩阵对于准确地估计机器人的状态和地图至关重要。
  • 矩阵的稀疏性表明大多数状态变量不直接相互影响,这在数值计算上是有利的。稀疏矩阵技术可以显著提高求解SLAM问题的效率。

因此,( A ) 矩阵不仅在物理上描绘了机器人在环境中的动态和相互作用,而且它是数学上解决状态估计问题的基础。通过理解和利用这种结构,可以更有效地实现机器人的自主导航和环境感知。

在SLAM(Simultaneous Localization and Mapping)问题中,矩阵 ( A ) 和向量 ( b ) 的物理意义及其在优化问题中的作用是核心部分,反映了状态估计的数学模型与实际操作间的联系。以下是使用 LaTeX 语法格式输出的详细解释:

矩阵 ( A ) 和向量 ( b ) 的物理意义

  1. 矩阵 ( A ):矩阵 ( A ) 在SLAM中代表系统的线性化状态转移和观测模型。它是状态和观测之间关系的线性近似表示,包含了如下两个主要部分:
    • 状态转移模型:通过 ( F_i^{i-1} ) 和 ( G_i^i ) 描述的块,这些块说明了给定控制输入 ( u_i ) 下,机器人状态从 ( x_{i-1} ) 转移到 ( x_i ) 的预测模型。
    • 观测模型:通过 ( H_k^{i_k} ) 和 ( J_k^{j_k} ) 描述的块,这些块表明如何从当前的状态和地标位置预测可能的观测结果 ( z_k )。
  2. 向量 ( b ):向量 ( b ) 通常包含观测数据与由模型预测的状态转移或观测的差异。在SLAM中,( b ) 可能由实际观测值 ( z_k ) 减去根据当前状态估计和地标估计得到的预测观测值构成。

最小化正规方程 ( A^T A x = A^T b )

在SLAM问题的线性化模型中,我们通常希望找到一个状态变量的更新量 ( ),使得模型预测尽可能贴近实际的观测。具体来说,我们尝试最小化以下目标函数:

\[ \| A \delta - b \|_2^2 \]

这里:

  • ( | |_2^2 ) 表示二范数的平方,用于度量误差的大小。
  • ( A ) 表示对当前状态估计的改正(或更新)。
  • ( b ) 是观测数据与模型预测之间的差异。

通过最小化这个误差,我们在数学上是在寻找一个向量 ( ),使得 ( A ) 尽可能接近 ( b ),即尝试对预测的状态和观测进行校正,以逼近真实的观测值。

解决这个最小化问题通常涉及到求解正规方程:

\[ A^T A \delta = A^T b \]

  • ( A^T A ) 是一个方阵,通常更易于处理,尤其是当 ( A ) 是稀疏或结构化的时,这种方法效率更高。
  • ( A^T b ) 是一个向量,表示在当前误差度量下,最佳更新方向的投影。

简而言之,矩阵 ( A ) 和向量 ( b ) 在SLAM中承载着将复杂的非线性系统模型化并使之适合数学优化的重任。通过这种优化,我们能够不断提高对机器人在未知环境中位置及其环境地图的估计精度。

因子图的迭代优化算法

图优化的迭代优化通常涉及到通过逐步优化来调整和改进图中的节点估计(即状态变量),以最小化整个图中所有误差函数(因子)的总和。这个过程通常使用类似于高斯-牛顿或者Levenberg-Marquardt(LM)算法的优化方法。

图优化迭代过程

  1. 初始化
    • 节点的初始估计通常基于传感器的原始数据或某些先验知识。
  2. 构建误差函数
    • 对于图中的每一个因子,构建一个误差函数,该函数衡量当前节点状态估计与观测数据之间的差异。
  3. 线性化误差函数
    • 在当前估计点附近对误差函数进行线性化,通常是通过计算误差函数相对于每个节点状态的雅可比矩阵。
  4. 构建正规方程
    • 将所有线性化的误差函数组合成一个大的系统方程。这通常表现为一个大的稀疏线性系统,可以写成 \[ H\delta = b \] 的形式,其中 ( H ) 是信息矩阵(或海森矩阵),而 ( b ) 是基于误差的向量。
  5. 求解线性系统
    • 使用直接求解器(如Cholesky分解)或迭代求解器(如共轭梯度法)来解决上述线性系统,以找到状态更新 ()。
  6. 更新状态估计
    • 将解得的 () 应用到当前的状态估计中,更新所有节点的状态。
  7. 迭代和收敛检查
    • 检查解的改进程度,如果未达到预设的收敛标准(如状态改变量的大小低于某个阈值或达到最大迭代次数),则重复步骤 2-6。
  8. 最终优化结果
    • 当迭代收敛时,当前的状态估计被认为是优化的解,即图中所有节点的最优估计。

使用Levenberg-Marquardt方法的特别说明

Levenberg-Marquardt算法是高斯-牛顿算法的一个变种,通过引入一个调整参数(阻尼因子)来平衡算法的迭代速度与稳定性。这个阻尼因子可以动态调整,以在必要时约束更新步长,从而提高收敛性和算法的鲁棒性。

图优化的迭代优化是一个强大的工具,特别适用于大规模、复杂的环境建模,如同时定位与地图创建(SLAM)等应用场景中。通过迭代改进,算法能够有效地精炼传感器数据与模型之间的不一致,最终达到高精度的状态估计。

为什么用Cholesky 分解求解正规方程

当然可以,这里是使用 LaTeX 语法重新整理的解释:

为什么使用 Cholesky 分解?

  1. 效率和稳定性
    • 效率:Cholesky分解相比于其他方法(如LU分解或QR分解)通常更为高效。它的计算复杂度大约是 \[\frac{1}{3} n^3\] 对于 \[n \times n\] 矩阵,而LU分解和QR分解的计算复杂度分别大约是 \[\frac{2}{3} n^3\]\[2n^3\]。这使得Cholesky分解成为求解大型正定线性系统的首选方法。
    • 稳定性:对于对称正定矩阵,Cholesky分解非常稳定。由于对称性和正定性保证了分解的成功执行,它避免了在数值分解过程中的潜在不稳定性。
  2. 数值属性
    • 在数值线性代数中,Cholesky分解对于对称正定矩阵来说是一种数值上稳定的方法,特别是当矩阵 \[A^T A\] 的条件数相对较低时。这是因为它直接利用了矩阵的对称性和正定性。
  3. 简化求解过程
    • 通过Cholesky分解,矩阵 \[A^T A\] 被分解为一个下三角矩阵 \[L\] 和其转置的上三角矩阵 \[L^T\] 的乘积,即 \[A^T A = L L^T\]。这种分解简化了线性系统的求解,因为只需先解下三角系统 \[L \mathbf{y} = A^T \mathbf{b}\] 获得临时向量 \[\mathbf{y}\],然后解上三角系统 \[L^T \mathbf{x} = \mathbf{y}\] 获取最终解 \[\mathbf{x}\]
    • 这种分步求解的过程比直接求解原方程更易于处理,并且可以有效地利用现代计算机体系结构。

总结

由于上述优点,Cholesky分解成为解决正定线性系统(特别是在涉及正规方程时)的理想选择,提供了一种高效、稳定的方法来求解最小二乘问题中的线性系统。这在实际应用中,如图优化、信号处理、统计数据分析等领域,具有广泛的应用。

为什么要引进\(G_i^i=-I_{d\times d}\)

在因子图和图优化中进行线性化处理时,引入矩阵\(G_i^i = -I_{d \times d}\)是为了便于表达和操作线性化的状态变化量。这里,$ I_{d d} $表示 ( d )-维的单位矩阵。矩阵 ( G_i^i ) 的作用是在数学表达式中简化状态变量 ( x_i ) 的运算,从而使得该线性系统易于使用线性代数技术求解。

背景

在线性化的SLAM或图优化问题中,我们通常处理的是两类项:

  1. 过程模型项:它们涉及从一个状态 ( x_{i-1} ) 到下一个状态 ( x_i ) 的过渡,如机器人的运动模型。
  2. 测量模型项:涉及从状态变量到测量变量的映射,如从机器人的位置到传感器观测的映射。

引入 $G_i^i = -I_{d d} $

在方程 (6) 中,我们考虑了状态变量 ( x_i ) 和 ( x_{i-1} ) 之间的关系。线性化后的模型需要处理这些状态变量的增量 ( x_i ) 和 ( x_{i-1} )。线性化过程后的表达式为:

\[ f_i(x_{i-1}, u_i) - x_i \approx \{ f_i(x_{i-1}^0, u_i) + F_i^{i-1} \delta x_{i-1} \} - \{ x_i^0 + \delta x_i \} \] 要将这个方程转换成标准的线性形式,我们需要确保 ( x_i ) 和 ( x_{i-1} ) 在同一方程中有适当的正负号。在这里,引入 ( G_i^i = -I_{d d} ) 允许我们将 ( x_i ) 项简单地表示为加法形式的负增量,即:

\[ F_i^{i-1} \delta x_{i-1} - \delta x_i = F_i^{i-1} \delta x_{i-1} + G_i^i \delta x_i \] 这样,线性方程更容易处理,因为它将状态变化量 ( x_i ) 直接以线性代数的标准形式表示,允许使用矩阵运算直接进行计算。

线性最小二乘问题

通过引入 ( G_i^i = -I_{d d} ),整个问题可表示为一个标准的线性最小二乘问题:

\[ \delta^* = \underset{\delta}{\operatorname*{argmin}} \left\{ \sum_{i=1}^M \| F_i^{i-1} \delta x_{i-1} + G_i^i \delta x_i + a_i \|_{\Lambda_i}^2 + \sum_{k=1}^K \| H_k^{i_k} \delta x_{i_k} + J_k^{j_k} \delta l_{j_k} - c_k \|_{\Sigma_k}^2 \right\} \] 这个形式使得应用线性代数和优化方法成为可能,从而可以有效求解大规模SLAM问题。

因子图的即插即用

因子图对于处理传感器更新频率不一致问题提供了一种灵活且有效的方式,特别是在复杂的多传感器融合系统中。这主要得益于因子图在处理多源信息与维持全局状态估计方面的灵活性。

在因子图中,每个传感器的输入可以被视为一个因子,它只与图中的某些变量(例如机器人的某个位置)相关联。当新的传感器数据到来时,可以简单地将相应的新因子添加到图中,而无需重新配置整个系统。这使得因子图非常适合于处理具有不同更新频率的传感器数据。

具体例子

假设一个机器人装备了一个高频率的IMU(惯性测量单元)和一个低频率的GPS。IMU能够提供高频率但低精度的位置和方向更新,而GPS提供低频率但高精度的位置测量。

  1. 构建因子图:机器人的每个状态(位置和姿态)都是图中的一个节点。IMU和GPS的每次读数分别产生一个新的因子,这些因子与特定的状态节点相连接。
  2. IMU因子:IMU的数据频繁更新,为图中相邻的状态节点间提供动态的约束(如速度和方向)。
  3. GPS因子:GPS数据不频繁,但每当数据到来时,会在因子图中创建一个新的因子,这个因子直接将GPS测量的位置与最近的状态节点相关联。
  4. 因子更新:无论何时收到新的传感器数据,都只需将对应的新因子添加到图中并重新优化,而不必修改或重处理整个图。

与卡尔曼滤波的比较

卡尔曼滤波器(包括其扩展和无迹变体)通常需要处理所有传感器数据的固定结构和预设的处理流程。每个传感器都需要在每个时间步被处理,即使某些传感器的数据并未更新。

  • 即插即用:卡尔曼滤波处理新传感器集成通常需要重新设计滤波过程和状态协方差的计算,而因子图可以更容易地通过添加新的因子来整合新传感器,无需重新设计或重新初始化整个系统。
  • 处理非同步数据:卡尔曼滤波器对于非同步传感器数据处理通常需要额外的插值或同步机制,而因子图可以更自然地通过在适当的时间点插入因子来处理这种数据。

因子图的这种灵活性使其非常适合于动态环境中的实时应用,特别是当涉及多个传感器,且这些传感器可能动态变化(如增加或移除传感器)时。因子图提供了一种高效、灵活且可扩展的方法来处理这种复杂性。

卡尔曼滤波器(KF)是一种适用于线性系统(或近似线性)的递推滤波方法,它依赖于数据输入的同步和连续性。在实际应用中,如自动驾驶车辆、机器人导航或多传感器融合系统中,传感器可能无法以同一频率或同步的方式产生数据。这会给基于卡尔曼滤波的系统设计带来挑战。

卡尔曼滤波器和非同步数据处理

当使用卡尔曼滤波器处理来自不同传感器的数据时,若这些传感器的更新频率不一致或数据非同步接收,可能需要进行额外处理以保证数据的时序一致性。这主要是因为卡尔曼滤波器在每个时间步更新状态估计时都假定接收到的数据是当前时刻的准确反映。

举例说明

假设有一个移动机器人使用两种传感器:一个激光雷达(LIDAR)和一个轮速计。激光雷达每100毫秒更新一次位置信息,而轮速计每30毫秒更新一次速度信息。如果直接将这些数据输入到一个卡尔曼滤波器,可能会出现问题,因为两种数据的时间戳不匹配。

为了同步这些传感器数据,一个常见的方法是使用插值:

  • 插值:例如,可以使用线性插值或更复杂的插值方法(如样条插值)来估计激光雷达测量间隙中的位置数据,或者在轮速计测量间隙中估计速度,以便在非测量时刻提供数据估计。
  • 数据缓冲:另一种方法是对快速更新的传感器(轮速计)数据进行缓冲,直到慢速传感器(激光雷达)的下一次数据到来,然后一起处理这些数据。

卡尔曼滤波器(KF)处理的核心在于它按照严格的时序顺序处理观测数据来更新状态估计。这个处理过程基于两个关键假设:系统的动态是已知的,并且观测数据是同步到特定的时间点的。当观测数据来自频率不同或非同步的多个传感器时,这些假设可能不再成立,直接接受数据而不进行任何处理可能会导致以下问题:

1. 误差累积

卡尔曼滤波器在每个时间步更新其状态估计时,假设最新的观测数据反映了当前的系统状态。如果直接使用非同步的观测数据,滤波器可能无法准确地关联观测数据与其对应的状态,从而导致误差累积。例如,如果一个传感器的数据延迟或提前了,而滤波器误以为这些数据是最新的,那么更新的状态估计可能会基于错误的信息。

卡尔曼滤波器(KF)处理的核心在于它按照严格的时序顺序处理观测数据来更新状态估计。这个处理过程基于两个关键假设:系统的动态是已知的,并且观测数据是同步到特定的时间点的。当观测数据来自频率不同或非同步的多个传感器时,这些假设可能不再成立,直接接受数据而不进行任何处理可能会导致以下问题:

1. 误差累积

卡尔曼滤波器在每个时间步更新其状态估计时,假设最新的观测数据反映了当前的系统状态。如果直接使用非同步的观测数据,滤波器可能无法准确地关联观测数据与其对应的状态,从而导致误差累积。例如,如果一个传感器的数据延迟或提前了,而滤波器误以为这些数据是最新的,那么更新的状态估计可能会基于错误的信息。

2. 时间戳不匹配

在处理来自不同传感器的数据时,每个传感器的输出通常带有时间戳。如果这些时间戳不一致,直接接受数据意味着滤波器可能会在错误的时间点应用观测数据,导致状态更新不准确。这对于那些高度依赖精确时间信息的应用(如导航和跟踪系统)尤为重要。

3. 数据丢失或重叠

非同步数据可能导致数据在时间轴上出现重叠或缺口。例如,一个传感器可能在另一个传感器的两次读数之间多次更新,或者在长时间内没有更新。这种情况下,直接处理可能导致某些观测数据被忽略或重复计算,影响滤波器的性能。

解决方法:插值和数据缓冲

为了解决这些问题,卡尔曼滤波器通常采用以下策略:

  • 插值:通过插值方法(如线性插值、样条插值等)创建虚拟的观测数据,填补时间上的缺口,使得每个时间步都有相应的数据输入。这有助于保持滤波器的连续性和准确性。

  • 数据缓冲:将来自不同传感器的数据存储在缓冲区中,直到收集到所有相关传感器的数据为止。然后,根据这些数据的时间戳一起处理它们。这样可以保证数据在时间上的一致性,减少由于时间错配引起的误差。

这两种方法使得卡尔曼滤波器能够更加稳定和准确地处理来自多源的非同步数据,从而提高系统的整体性能。

在处理来自不同传感器的数据时,每个传感器的输出通常带有时间戳。如果这些时间戳不一致,直接接受数据意味着滤波器可能会在错误的时间点应用观测数据,导致状态更新不准确。这对于那些高度依赖精确时间信息的应用(如导航和跟踪系统)尤为重要。

3. 数据丢失或重叠

非同步数据可能导致数据在时间轴上出现重叠或缺口。例如,一个传感器可能在另一个传感器的两次读数之间多次更新,或者在长时间内没有更新。这种情况下,直接处理可能导致某些观测数据被忽略或重复计算,影响滤波器的性能。

解决方法:插值和数据缓冲

为了解决这些问题,卡尔曼滤波器通常采用以下策略:

  • 插值:通过插值方法(如线性插值、样条插值等)创建虚拟的观测数据,填补时间上的缺口,使得每个时间步都有相应的数据输入。这有助于保持滤波器的连续性和准确性。

  • 数据缓冲:将来自不同传感器的数据存储在缓冲区中,直到收集到所有相关传感器的数据为止。然后,根据这些数据的时间戳一起处理它们。这样可以保证数据在时间上的一致性,减少由于时间错配引起的误差。

这两种方法使得卡尔曼滤波器能够更加稳定和准确地处理来自多源的非同步数据,从而提高系统的整体性能。

因子图的优势

相比之下,因子图提供了一个更灵活的框架来处理非同步数据。因子图是一种图形模型,它以图的形式表示变量之间的条件依赖关系。在因子图中,每个传感器的读数可以作为一个因子被添加到图中,这些因子与它们影响的状态变量相关联。如果传感器数据是非同步的,可以:

  • 动态添加因子:根据每个传感器的时间戳动态地向图中添加新的因子。这意味着每个传感器更新只影响相关的状态变量,不需要人为地同步所有数据。
  • 时间标记:每个因子可以携带时间信息,因此状态更新可以在正确的时间点进行,与数据实际测量的时间相匹配。

因此,因子图模型通过在合适的时间点添加因子来自然地处理非同步的传感器数据,从而避免了在卡尔曼滤波器中必须进行的数据插值或同步。这样不仅简化了处理过程,而且可能提高了系统的整体性能和精度。

因子图优化更新节点状态

  • 应用方式:图优化通常处理的是整个数据集(如历史路径或多个观测数据点),并且通常适用于离线或批处理模式,适用于大规模空间问题,如SLAM。而EKF通常在线性化当前或最近的状态来进行实时更新,适用于需要快速反应的动态系统,如实时机器人导航。
  • 数据处理方式:EKF以递归方式连续更新状态,侧重于实时估计;图优化则可能在获得新数据后对整个历史数据进行重新优化,侧重于精度和全局一致性。

让我们逐一探讨这两个问题:

1. 根据导数信息计算状态更新量的例子

考虑一个简单的SLAM问题,其中机器人需要确定它的位置 ( x ) 和地图上一个地标的位置 ( l )。假设有一个观测 ( z ) 与机器人的位置和地标之间的距离相关。

步骤

  • 误差函数:设 ( h(x, l) ) 是从机器人位置 ( x ) 到地标 ( l ) 的预测测量,误差函数 ( e ) 可表示为 ( e = z - h(x, l) )。

  • 雅可比矩阵:计算误差 ( e ) 关于状态 ( x ) 和地标位置 ( l ) 的雅可比矩阵 ( J_x ) 和 ( J_l )。

  • 更新公式:使用高斯牛顿法,更新公式可以写作:

    \[ \begin{pmatrix} \Delta x \\ \Delta l \end{pmatrix} = -(J^T J + \lambda I)^{-1} J^T e \]

    其中 ( J ) 是组合了 ( J_x ) 和 ( J_l ) 的雅可比矩阵,( ) 是为了提高数值稳定性而添加的阻尼因子(Levenberg-Marquardt算法中使用),( I ) 是单位矩阵。

  • 状态更新:计算得到的 ( x ) 和 ( l ) 用来更新当前状态和地标位置。

2. 增量式图优化方法的实际操作

在实时系统如自动驾驶车辆中,增量图优化主要关注最近的观测和状态,以减少每次更新所需的计算资源。

步骤

  • 局部图更新:每当接收到新的观测数据时,只更新与这些新数据直接相关的节点和边,这些通常是最近的机器人状态和相关地标。
  • 窗口优化:一种常见的策略是固定滑动窗口方法,即仅在固定大小的最近状态窗口内进行优化。较旧的状态和地标,一旦它们超出窗口,可以被边缘化(从优化问题中移除但考虑其对现有状态的统计影响)或锁定。
  • 增量求解器:使用专门设计的增量求解器,如 iSAM 或 g2o,这些工具能够有效地更新信息矩阵和求解增量更新问题。
  • 实时反馈:通过快速更新仅涉及局部改动的解决方案,系统能够快速地对新数据做出响应。

这种增量更新方法能够显著降低因每个新观测而进行全局优化的计算成本,使系统更适合实时应用场景,如自动驾驶或动态环境中的机器人导航。

如何根据联合概率密度分布 求解

在解决导航和地图构建的问题中,目标通常是找到一组变量(机器人的位置和地图配置)的配置,这些配置最能解释观测数据和已知信息。这里的“解释”通常是指在给定模型下,找到最大化联合概率 \(P(X, L, Z)\) 的 (X)(位置)和 (L)(地图)的配置。这种方法是基于最大似然估计或最大后验估计的原理。

最大似然估计和最大后验估计

  • 最大似然估计(MLE): 在这种情况下,我们寻找的是最大化观测数据 (Z) 出现概率的位置 (X) 和地图 (L) 配置,不考虑位置和地图的先验概率。数学上表达为 \(\arg\max_{X, L} P(Z|X, L)\)

  • 最大后验估计(MAP): 更常见的是,我们寻找能够最大化后验概率\(P(X, L | Z)\)的 (X) 和 (L)。由于根据贝叶斯公式,\(P(X, L | Z) \propto P(Z | X, L) P(X, L)\),这不仅考虑了观测数据的似然度,还包括了关于位置和地图的先验知识。

计算过程

  1. 定义概率模型:
    • 状态转移概率 \(P(x_i | x_{i-1}, u_i)\)描述了给定前一状态和控制输入后,当前状态的概率。
    • 观测模型 \(P(z_k | x_{i_k}, l_{j_k})\) 描述了给定当前状态和地图配置时,观测数据出现的概率。
    • 先验\(P(X, L)\)可能基于以前的知识或假设(如位置的均匀分布,或地图的某些特征)。
  2. 优化过程:
    • 在实践中,最大化\(P(X, L | Z)\) 通常涉及到优化算法,如梯度上升、期望最大化(EM)算法或其他非线性优化技术。
    • 在机器人技术中,常见的优化工具包括粒子滤波(用于非线性/非高斯问题)和卡尔曼滤波(用于线性/高斯问题)。
  3. 联合概率和解释:
    • 我们不是寻找最小的联合概率,而是寻找最大化后验概率的 (X, L) 配置,因为后验概率高意味着给定观测数据和先验知识下,该配置更可能是正确的。

结论

通过最大化后验概率\(P(X, L | Z)\),我们找到一组 (X) 和 (L),使得在已有的观测 (Z) 和可能的先验知识下,这组配置是最有可能的。这就是如何使用联合概率模型来优化并解释机器人的导航和地图构建问题。这种方法使得可以在不确定的环境中做出最好的估计和决策。

为什么用联合概率密度来描述导航问题

在导航和其他机器人感知系统中,状态估计问题经常涉及位置(X),地图(L)和观测数据(Z)的不确定性。利用联合概率密度函数来表示这些变量,可以系统地包含所有相关的不确定性和依赖关系,从而使得可以进行精确的推断和决策。以下是详细解释为什么采用这种表示方法:

在导航问题中:

  • X 表示机器人的位置状态序列,即在不同时间点的位置状态。
  • L 表示环境的地图,包括可能的地标或其他特征。
  • Z 表示机器人的观测数据序列,即机器人通过其传感器相对于环境特征的观测。

概率模型表示

  1. 初始状态的概率 (P(x_0)):这是机器人在开始导航前的位置状态的初始概率分布。这通常是基于先验知识或初始观测所给定的。

  2. 状态转移概率 ({i=1}^M P(x_i|x{i-1}, u_i)):这一部分表示机器人位置的演化,其中每个状态 (x_i) 取决于前一个状态 (x_{i-1}) 和执行的控制指令 (u_i)。状态转移概率捕捉了系统动态和可能的过程噪声。

  3. 观测模型 ({k=1}^K P(z_k|x{i_k}, l_{j_k})):每个观测 (z_k) 都依赖于特定时刻 (i_k) 的位置 (x_{i_k}) 和相关联的地图特征 (l_{j_k})。观测模型描述了如何根据机器人的当前位置和周围环境的地标特征产生观测数据,包括观测噪声。

联合概率密度函数的优势

使用联合概率密度函数 (P(X, L, Z)) 的优势在于:

  • 全面性:通过整合所有可能的不确定性来源(控制输入的不确定性、位置的不确定性和观测的不确定性),这种方法提供了一个全面的系统描述。
  • 条件依赖性的显式表达:明确表达了各状态和观测之间的条件依赖关系,这对于理解和实现滤波器和估计算法至关重要。
  • 灵活性:这种表示允许使用各种概率推断技术,如卡尔曼滤波、粒子滤波或图优化,来估计位置和地图。
  • 信息整合:可以有效地整合多时刻和多来源的信息,优化状态估计和地图构建。

总之,联合概率密度函数通过对所有变量的依赖关系进行建模,为解决导航中的状态估计问题提供了一个强大的数学框架,使得可以利用概率方法处理不确定性和推导最优解。这种方法是解决SLAM问题的基石,也是现代自动驾驶和机器人导航系统中不可或缺的一部分。

再尝试用一个更简单和直观的方式来理解为什么在导航中的状态估计问题可以通过一个联合概率密度函数来表示。我们可以通过具体的例子和逐步解析来帮助理解。

一个具体的例子

假设我们有一个机器人在一个未知的环境中导航。机器人的任务是理解自己的位置(状态X),探索并绘制周围的地图(L),同时根据自己的传感器收集信息(观测Z)。

  1. 初始状态 (P(x_0)):
    • 在一开始,我们只知道机器人可能在一个初始位置,比如我们可以假设它以某个概率分布(如高斯分布)位于起点附近。
  2. 状态转移 (P(x_i|x_{i-1}, u_i)):
    • 每次机器人移动,它的新位置 (x_i) 将基于前一个位置 (x_{i-1}) 和其执行的移动指令 (u_i)(比如向前移动一米)。这个过程不是完全准确的,因此我们使用概率模型(如另一个高斯分布)来表示可能的新位置。
  3. 观测模型 (P(z_k|x_{i_k}, l_{j_k})):
    • 同时,机器人会使用传感器(如摄像头或激光雷达)观测周围环境。这些观测 (z_k) 取决于机器人的当前位置 (x_{i_k}) 和它所观测到的环境特征 (l_{j_k})(比如附近的一个门或窗户)。观测也不是完全精确的,同样用概率模型来描述。

为什么使用联合概率密度函数

当我们想要同时考虑所有这些因素来估计机器人的最可能的位置、绘制地图和解释观测数据时,我们可以构造一个联合概率模型 (P(X, L, Z)):

  • 综合所有信息:这个联合模型把所有关于机器人位置的不确定性、地图的不确定性和观测数据的不确定性都整合在一起。我们通过计算这个联合概率来寻找最有可能解释所有观测数据和已知信息的位置和地图配置。
  • 条件依赖关系:这个模型同时还能清晰地表达各种条件依赖关系,例如观测依赖于位置和地图的某些特定特征,位置的更新依赖于先前的位置和移动指令。
  • 灵活的推断和优化:有了这个联合模型,我们可以利用不同的数学工具和算法(如贝叶斯滤波器、粒子滤波、最大似然估计等)来推断未知的状态和地图。

总结来说,这个联合概率密度函数提供了一种全面并且系统的方法来考虑和处理导航问题中的所有不确定性和依赖性。这种方法使得可以更有效地进行状态估计和环境映射,帮助机器人更好地理解和操作在其所处的环境中。

因子图通过将测量值作为已知的参数而不是变量来简化问题表示?

因子图通过将测量值作为已知的参数而不是变量来简化问题表示,使得模型更直接地关联到优化问题本身,这对于高效解决SLAM等复杂问题至关重要。 这段话怎么理解?通过将测量值作为已知的参数而不是变量来简化问题表示?

这段话指出因子图在处理如SLAM这样的问题时的一个关键特性:将测量数据视为已知参数,而不是作为需要求解的变量。这种处理方式在建模和求解复杂系统时有几个重要的好处:

简化问题表示

  1. 明确角色分配:在SLAM等问题中,我们的主要目标是推断或优化未知的状态变量(如机器人的位置、地图的特征点等),而测量数据(如来自传感器的读数)是用来帮助确定这些未知变量的。将测量值视为已知参数,意味着这些数据在模型中的角色是辅助性的,它们用于构建关于未知变量的条件概率表达,而不是自身需要被解决的问题。

  2. 减少求解复杂度:如果将测量值也视为变量,那么问题的复杂度会显著增加,因为这样不仅需要推断未知的状态变量,还要同时解决这些测量值的准确性和可靠性问题。将测量值作为参数,可以直接使用这些数据来形成对未知变量的概率描述,从而简化了模型和计算。

直接关联到优化问题本身

  1. 更加聚焦核心任务:在SLAM等应用中,核心任务是基于测量数据推断或更新未知的状态变量。将测量数据视为参数,可以使模型更加集中于如何利用这些数据来改善对未知变量的估计,而非分散注意力于处理测量数据本身的不确定性。

  2. 提高算法效率:当模型聚焦于核心变量时,优化算法可以更有效地调整和优化这些变量的值以达到最佳估计。这通常会导致更快的计算速度和更好的优化性能,因为算法可以直接利用固定的、已知的测量数据来逐步改进状态估计。

总之,将测量值作为已知参数处理,有助于因子图更加高效和直接地应对复杂的优化问题,如SLAM,这样做能够减轻计算负担,同时提高问题求解的准确性和效率。这种方法使得模型设计和算法实现都更加直接和实用,特别是在实时或资源受限的应用场景中。