用于误差状态卡尔曼滤波的四元数运动学
摘要
本文全面修订了与四元数及三维空间旋转相关的概念和公式,并详述了它们在误差状态卡尔曼滤波器等估计引擎中的正确使用方法。
该论文深入研究了旋转群及其李群结构,同时采用四元数和旋转矩阵两种表述形式。特别关注了旋转扰动、导数与积分的定义,并提供了大量直观解释和几何诠释,以帮助读者理解三维旋转的内在机制。
本材料旨在设计适用于实际应用的误差状态卡尔曼滤波器的精确公式,通过整合来自惯性测量单元(IMU)的信号实现。
1 四元数的定义与性质
1.1 四元数的定义
关于四元数的一个特别吸引我的介绍方式是通过凯莱-迪克森构造给出的:
若有两个复数 和 ,通过构造 并定义 ,便可得到一个四元数空间 中的数,
| (1) |
其中 和 是三个虚数单位,其定义满足
| (2a) | |||
| from which we can derive | |||
| (2b) | |||
由(1)式可见,我们可以在四元数定义中嵌入复数,进而嵌入实数和纯虚数,这意味着实数、虚数和复数本质上都是四元数,
| (3) |
同样地,为了完备性起见,我们可以在 的三维虚子空间中定义数。我们将其称为纯四元数 ,并可用 表示纯四元数空间,
| (4) |
值得注意的是,虽然单位长度的常规复数 可通过单一复数乘法 编码二维平面中的旋转,但"扩展复数"或单位长度的四元数 则能通过双四元数乘法 (如本文后续所述)编码三维空间中的旋转。
注意: 并非所有四元数定义都相同。部分作者将乘积写作 而非 ,因此他们得到性质 ,这导致 以及左手系四元数。此外,许多作者将实部置于末位,形成 。这些选择并无本质区别,但会使整个公式体系在细节上呈现差异。更多解释与歧义消除请参阅第 3 节。
注意: 还存在其他使公式细节产生差异的约定。这些约定涉及我们对旋转算子"含义"或"解释"的理解——无论是旋转矢量还是旋转参考系(本质上构成互逆操作)。更多解释与歧义消除请同样参阅第 3 节。
1.1.1 四元数的替代表示法
实数+虚数的表示法 并不总是适合我们的需求。在运用代数(2)的前提下,四元数可表示为标量+向量的组合,
| (5) |
其中 被称为实部或标量部 , 则称为虚部或向量部 。 111Our choice for the subscripts notation comes from the fact that we are interested in the geometric properties of the quaternion in the 3D Cartesian space.
Other texts often use alternative subscripts such as or , perhaps better suited for mathematical interpretations.
亦可将其定义为标量-向量的有序对
| (6) |
我们通常将四元数 表示为四维向量 ,
| (7) |
这使得我们能够使用矩阵代数来处理涉及四元数的运算。
在某些情况下,我们可以通过滥用“ ”符号来混合使用不同的表示法。典型的例子包括实四元数和纯四元数 ,
| (8) |
1.2 四元数主要性质
1.2.1 加法
该求和过程直接明了,
| (9) |
根据构造,该和具有交换律与结合律性质,
| (10) | ||||
| (11) |
1.2.2 四元数乘积
记作 的四元数乘积运算需使用原始形式(1)及四元数代数规则(2),将其写成向量形式可得
| (12) |
该运算亦可表述为标量部分与向量部分的组合:
| (13) |
其中叉积的存在表明,四元数乘法在一般情况下是不可交换的 ,
| (14) |
这种普遍非交换性的例外情况仅限于以下情形:当其中一个四元数为实数 时,或 ,或 时;或者当两个向量部分平行时 。只有在这些情况下,四元数乘法才是可交换的。
然而,四元数乘法具有结合律特性,
| (15) |
并且对和具有分配性 ,
| (16) |
两个四元数的乘积是双线性的,可以表示为两种等价的矩阵乘积形式,即
| (17) |
其中 和 分别为左、右四元数乘积矩阵,通过简单观察即可由(12)和(17)式导出,
| (18) |
| (19) |
此处, 斜算子 222The skew-operator can be found in the literature in a number of different names and notations, either related to the cross operator , or to the ‘hat’ operator ∧, so that all the forms below are equivalent,
生成叉积矩阵,
| (20) |
该矩阵为斜对称矩阵 ,等价于叉积运算, 即 ,
| (21) |
配备乘积运算 的四元数构成一个非交换群。该群的单位元 和逆元 将在下文探讨。
1.2.3 单位元
关于该乘积运算的单位四元数 满足 。它对应于实数乘积中的单位元"1",以四元数形式表示,
1.2.4 共轭
四元数的共轭定义为
| (23) |
其具有以下性质
| (24) |
且
| (25) |
1.2.5 范数
四元数的范数定义为
| (26) |
它具有以下性质
| (27) |
•
1.2.6 逆
四元数的逆 满足该四元数与其逆相乘等于单位元,
| (28) |
其计算方式为
| (29) |
1.2.7 单位或归一化四元数
对于单位四元数, ,因此
| (30) |
当将单位四元数解释为方向规范或旋转算子时,该性质意味着逆旋转可通过共轭四元数实现。单位四元数总可表示为以下形式:
| (31) |
其中 是单位向量, 是标量。
根据(27),赋予乘积运算 的单位四元数构成一个非交换群,其中逆元与共轭元重合。
1.3 四元数的附加性质
1.3.1 四元数交换子
1.3.2 纯四元数的乘积
1.3.3 纯四元数的自然幂
1.3.4 纯四元数的指数
四元数指数函数是与普通指数函数类似的四元数函数。与实数指数情形完全一致,其定义为绝对收敛的幂级数:
| (39) |
显然,实四元数的指数函数与普通指数函数完全一致。
更有趣的是,纯四元数 的指数定义为一个新的四元数,
| (40) |
令 ,其中 和 为单位量,并考虑(37),我们将级数中的标量项与向量项分组,
| (41) |
并分别识别出它们对应于 和 的级数展开。 333We remind that , and .
由此可得
| (42) |
这构成了欧拉公式的一个优美扩展, ,该公式最初是为虚数定义的。注意到由于 ,纯四元数的指数是一个单位四元数。同时还需注意这一性质,
| (43) |
对于小角度四元数,我们通过展开 和 的泰勒级数并截断来避免 中的除以零问题,从而获得不同精度的近似解。
| (44) |
1.3.5 一般四元数的指数
1.3.6 单位四元数的对数
显然可见,若 ,
| (47) |
即,单位四元数的对数是一个纯四元数。通过反转(42)可以轻松获得角度-轴值,
| (48) | ||||
| (49) |
对于小角度四元数,我们通过表达 的泰勒级数并截断, 444We remind that , and . 获得不同阶数的近似值,以避免除以零的情况,
| (50) |
1.3.7 一般四元数的对数
进一步地,若 为一普通四元数,
| (51) |
1.3.8 类型 的指数形式
2 旋转与交叉关系
2.1 三维向量旋转公式
图 1: 向量 绕轴 旋转角度 的示意图。详见正文说明。
我们在图 1 中展示了遵循右手定则的三维向量 绕单位向量 所定义轴线旋转角度 的过程。这一旋转通过将向量 分解为平行于 的分量 和正交于 的分量 来实现,即
这些分量可简便计算得出(其中 表示向量 与轴线 之间的夹角),
旋转时,平行分量不发生转动,
而正交分量将在与 垂直的平面内发生旋转。具体而言,若我们构建该平面的正交基 ,满足
条件 ,则有 。在该平面上进行 弧度的旋转将产生,
其发展过程为,
加上平行部分后,便得到了旋转向量的表达式 ,这被称为向量旋转公式
| (54) |
2.2 旋转群
在文献[0]中,旋转群[1]定义为原点旋转在复合运算下构成的群。旋转是保持向量长度和相对向量方向(即手性)的线性变换。其在机器人学中的重要性在于它表示三维空间中刚体的旋转:刚体运动要求运动过程中刚体内部的距离、角度和相对方向必须保持不变——否则,若范数、角度或相对方向未能保持,该物体就不能被视为刚体。
让我们通过满足这些性质的算子来定义旋转。
作用于向量[1]的旋转算子[0]可以从欧几里得空间的度量(由点积和叉积构成)定义如下。
-
•
旋转变换保持向量范数不变,(55a) -
•
旋转保持向量间的角度不变(55b) -
•
旋转保持了向量间的相对方向关系,(56)
容易证明前两个条件是等价的。因此我们可以将旋转群[0]定义为,
| (57) |
旋转群通常由旋转矩阵的集合来表示。然而,四元数同样能很好地表示它。本章旨在证明这两种表示方式具有同等有效性。正如读者将在表 1 中所见,它们在概念和代数层面都展现出诸多相似性。
表 1: 用于表示 的旋转矩阵与四元数
| 旋转矩阵, | 四元数, | |
| 参数 | ||
| 自由度 | 3 | 3 |
| 约束条件 | ||
| 约束条件 | ||
| ODE | ||
| 指数映射 | ||
| 对数映射 | ||
| 与 的关系 | 单一覆盖 | 双层覆盖 |
| 恒等元 | ||
| 逆元 | ||
| 复合运算 | ||
| 旋转算子 | ||
| 旋转作用 | ||
| 插值 | ||
| 交叉关系 | ||
或许最重要的区别在于,单位四元数群构成了 的双重覆盖(因此严格来说并非 本身),这一特性在我们的大多数应用中并不关键。 555The effect of the double cover needs to be considered when performing interpolation in the space of rotations. This is however easy, as we will see in Section 2.7.
为便于快速比较和评估,表格已前置插入。
的旋转矩阵与四元数表示形式将在后续章节中探讨。
2.3 旋转群与旋转矩阵
由于 运算符是通过标量积和向量积(二者均为线性运算)定义的,因此它是线性的。该运算可通过矩阵 表示,并通过矩阵乘积对向量 实施旋转。
| (58) |
将其代入(55a),利用点积 并展开可得对于所有 ,
| (59) |
由此得到关于 的正交性条件,
| (60) |
上述条件实为正交性条件,因为通过令 并代入上式可观察到,矩阵 的列向量 (其中 )具有单位长度且彼此正交
保持向量范数与角度不变的变换集合因此被称为正交群 ,记作 。正交群包含旋转(刚体运动)与反射(非刚体运动)。此处群的概念本质上(非正式地)意味着两个正交矩阵的乘积仍是正交矩阵 666Let and be orthogonal,
and build .
Then .
且每个正交矩阵都存在逆矩阵。
实际上,正交性条件(60)意味着逆旋转可通过转置矩阵实现,
| (61) |
添加相对方向条件(56)可保证刚体运动(从而排除反射),并对 、 777Notice that reflections satisfy , and do not form a group since . 产生一个额外约束
| (62) |
具有正单位行列式的正交矩阵通常被称为真或特殊正交矩阵。这类特殊正交矩阵的集合构成 的一个子群,称为特殊正交群 。作为群结构,两个旋转矩阵的乘积始终是旋转矩阵。 888See footnote 6 for and add this for : let , then .
2.3.1 指数映射
指数映射(以及下节将介绍的对数映射)是在三维旋转空间中实现便捷严谨运算的强大数学工具。它开启了适用于旋转空间的微积分体系之门,使我们能正确定义导数、扰动与速度并进行运算。因此该映射在旋转/姿态空间的估计问题中具有基础性地位。
旋转构成刚体运动。这种刚性意味着可以在 中定义一条连续轨迹或路径 ,将刚体从其初始方位 连续旋转至当前方位 。由于连续性,研究此类变换的时间导数是合理的。我们通过推导刚见过的性质(60)和(62)来实现这一点。
首先注意到,在满足 (60)条件的同时,不可能连续脱离单位行列式条件 (62),因为这将导致行列式从 跳跃至 。 999Put otherwise: a rotation cannot become a reflection through a continuous transformation.
因此,我们只需研究正交性条件(60)的时间导数。其表达式为
| (63) |
由此可得
| (64) |
这意味着矩阵 是斜对称的( 即 ,等于其转置的负矩阵)。斜对称 矩阵的集合记为 ,称为 的李代数 。斜对称 矩阵具有如下形式:
| (65) |
它们具有 3 个自由度,对应于我们已在(20)中介绍的叉积矩阵。这建立了一一映射 。现取向量 并写作
| (66) |
由此导出常微分方程(ODE)
| (67) |
在原点附近,我们有 ,上述方程简化为 。因此,可将李代数 解释为 在原点处的导数空间;它构成了 的切空间 ,或称速度空间 。基于这些事实,我们完全可以将 称为瞬时角速度矢量。
若 为常数,上述微分方程可进行时间积分如下
| (68) |
其中指数 由其泰勒级数定义,如下节所示。由于 和 是旋转矩阵,显然 也是一个旋转矩阵。定义向量 为编码整个周期 内完整旋转的旋转向量,我们得到
| (69) |
这被称为指数映射,是从 到 的一个映射,
| (70) |
2.3.2 大写指数映射
上述指数映射有时会以符号滥用的方式表示, 即将 与 混淆。为避免可能的歧义,我们选择用大写 的显式符号表示这个新映射 ,如(见图 2)所示
图 2: 旋转矩阵的指数映射
| (71) |
它与指数映射的关系是显而易见的,
| (72) |
在接下来的章节中,我们将看到被称为旋转向量或角度-轴向量的向量 ,通过 编码了旋转角度 和旋转轴 。
2.3.3 旋转矩阵与旋转向量:罗德里格斯旋转公式
2.3.4 对数映射
我们定义对数映射为指数映射的逆运算,
| (79) |
随着
| (80) | ||||
| (81) |
其中 是 的逆运算,即满足 与 关系式。
我们还定义了其大写形式 ,通过该形式可直接从旋转矩阵反解出旋转矢量 ,
| (82a) | |||
它与对数映射的关系是简单的,
| (83) |
2.3.5 旋转作用
2.4 旋转群与四元数
出于教学目的,我们着重强调四元数与旋转矩阵作为旋转群 表示形式之间的关联性
为此,我们最初将著名的四元数旋转作用公式作为假设,该公式表述为:
| (86) |
这使得我们能够通过重述旋转矩阵所用论述的方式,完整地展开四元数章节。这一假设的准确性将在稍后的 2.4.5 节中得到证明,从而验证该方法的有效性。
让我们将上述旋转代入正交性条件(55a),并利用(27)展开得到
| (87) |
由此可得 ,即四元数的单位范数条件,其表达式为
| (88) |
该条件类似于我们在旋转矩阵中遇到的(60)条件,其表达式为 。建议读者稍作停留,思考二者的相似之处。
单位四元数集合在乘法运算下构成一个群。该群在拓扑结构上是三维球面,即四维单位球体的三维表面 ,通常记作 。
2.4.1 指数映射
考虑一个单位四元数 ,即满足 ,我们按照处理旋转矩阵正交条件 的方式进行推导。对其求时间导数,
| (90) |
由此可得
| (91) |
这意味着 是一个纯四元数( 即 ,它等于其共轭的相反数,因此实部为零)。因此我们取一个纯四元数 并表示为,
| (92) |
左乘 得到微分方程,
| (93) |
在原点附近,我们有 ,上述方程简化为 。因此,纯四元数的空间 构成了四元数单位球面 的切空间 ,即其李代数。然而在四元数情形下,这个空间并非直接对应于速度空间,而是半速度空间,这一点我们稍后将详细说明。
若 为常数,该微分方程可积分得到
| (94) |
由于 和 均为单位四元数,其指数映射 仍为单位四元数——这从四元数指数性质(42)中已可推知。定义 后可得
| (95) |
这仍是一个指数映射:从纯四元数空间到由单位四元数表征的旋转空间的映射
| (96) |
2.4.2 大写指数映射
如我们所见,指数映射(96)中的纯四元数 通过 编码了旋转轴 与旋转角度的一半 。关于这个半角特性,我们将在 2.4.5、2.4.6 和 2.8 节详细解释。目前只需说明:由于旋转作用是通过双乘积 实现的,向量 经历的旋转是 编码角度的"两倍",或者说四元数 编码了对 预期旋转的"一半"。
为了直接表达角轴旋转参数 与四元数之间的关系,我们定义了一个大写的指数映射版本,该映射捕捉了半角效应(见图 3),
图 3: 四元数的指数映射。
| (97) |
其与指数映射的关系是显而易见的,
| (98) |
2.4.3 四元数与旋转向量
2.4.4 对数映射
我们还定义了大写对数映射,它能直接在笛卡尔三维空间中提供旋转角度 和旋转轴 ,
| (103) |
它与对数映射的关系是简单的,
| (104) |
2.4.5 旋转作用
2.4.6 流形的双重覆盖
考虑单位四元数 。当视为常规4维向量时, 与表示方向原点的单位四元数 之间的夹角 为,
| (110) |
同时,四元数 在三维空间中对物体旋转的角度 满足
| (111) |
这意味着我们有 ,因此四元数向量与四维空间中恒等元之间的夹角,等于该四元数在三维空间中所旋转角度的一半,
| (112) |
我们在图 4 中展示了这种双重覆盖现象。当两个四元数向量之间的夹角达到 时,三维旋转已完成 (即半周旋转);而当四元数向量完成半周旋转 时,三维旋转已实现完整一周旋转。四元数向量的后半周旋转 对应着三维旋转的第二个完整周期 ,即旋转流形的第二次覆盖。
图 4: 旋转流形的双重覆盖。左:单位 3-球面上的四元数 与单位四元数 形成角度 。中:所得 3D 旋转 的角度 是原始四元数的两倍。右:将 4D 和 3D 旋转平面叠加后可见,四元数 在 3-球面上旋转一周(红色)对应被旋转向量 在 3D 空间中旋转两周(蓝色)。
2.5 旋转矩阵与四元数
正如我们所见,给定旋转向量 ,其指数映射
单位四元数
和
旋转矩阵
产生旋转算子
并且
它们以相同的角度 绕同一轴 旋转向量 。 101010The obvious notation ambiguity between the exponential maps and is easily resolved by the context:
at occasions it is just the type of the returned value, or ;
other times it is the presence or absence of the quaternion product .
也就是说,如果
| (113) |
那么,
| (114) |
由于该恒等式两边对 呈线性关系,通过展开左侧并匹配右侧项,可得到与四元数等价的旋转矩阵表达式,即四元数转旋转矩阵公式:
| (115) |
本文中以 表示。四元数乘积的矩阵形式(17–19)提供了另一种计算公式,因为
| (116) |
经过一些简单推导后可得
| (117) |
旋转矩阵 相对于四元数具有以下性质:
| (118) | ||||
| (119) | ||||
| (120) | ||||
| (121) |
其中我们观察到:
(118) 单位四元数表示零旋转;(119) 四元数与其负数表示相同旋转,构成 的双重覆盖;(120) 共轭四元数表示逆旋转;且(121) 四元数乘积按与旋转矩阵相同的顺序组合连续旋转。
此外,我们有以下性质
| (122) |
该关系描述了四元数与旋转矩阵在运行标量 上的球面插值。
2.6 旋转合成
图 5: 旋转组合。在 中我们只需执行 操作,其中"求和"运算是可交换的。而在 中组合运算满足 ,其矩阵形式为 。这些算子不可交换,必须严格遵循顺序——恰当的符号表示有助于理解:'AB'与'BC'串联形成'AC'。
关于符号表示的说明
恰当的符号表示有助于确定复合运算中因子的正确顺序,特别是对于多次旋转的复合运算(参见图 5)。例如,令 (相应地 )表示从状态 到状态 的旋转,即 (相应地 )。那么给定由四元数 表示的一系列旋转时,我们只需将索引串联起来即可得到:
且已知逆向旋转对应于共轭运算 或转置运算 ,因此还可得出:
2.7 球面线性插值(SLERP)
四元数非常适合用于计算精确的姿态插值。
给定由四元数 和 表示的两种姿态,我们希望找到一个四元数函数 ,能够实现从 到 的线性插值。这种插值的特点是:当 从 变化到 时,物体会以恒定角速度绕固定轴,从姿态 连续旋转至姿态 。
方法一
第一种方法采用四元数代数,并遵循文献 中的几何推理,这应与目前所述内容易于关联。首先,计算从 到 的姿态增量 ,使得 ,
| (124) |
随后利用对数映射 111111We can use here either the maps and , or their capitalized forms and . The involved factor 2 in the resulting angles is finally irrelevant as it cancels out in the final formula. 求得对应的旋转向量 ,
| (125) |
最后保持旋转轴 不变,取旋转角度 的线性分数部分。通过指数映射 将其转化为四元数形式,并与原始四元数组合得到插值结果,
| (126) |
整个过程可表示为 ,简化为
| (127) |
并且通常(参见 53)被实现为:
| (128) |
图 6: 在 单位球面上的四元数插值,以及 旋转平面 上的情况正视图。
注:
方法二
还可开发其他独立于四元数代数内部结构、甚至独立于弧线嵌入空间维度的 Slerp 方法。
特别参见图 6,我们可以将四元数 和 视为单位球面上的两个单位向量,并在同一空间中进行插值。插值结果 是以恒定角速度沿连接 到 的最短球面路径的单位向量。该路径是通过单位球面与由 、 及原点定义的平面(图中虚线圆周)相交形成的平面弧线。关于这些方法与前述等价性的证明,参见 Dam 等人(1998)。
第一种方法采用向量代数,直接遵循上述思路。
将 和 视为两个单位向量;它们之间的夹角 121212The angle is the angle between the two quaternion vectors in Euclidean 4-space, not the real rotated angle in 3D space, which from (125) is . See Section 2.4.6 for further details.
可通过标量积导出,
| (131) |
我们按如下步骤进行。
我们确定旋转平面(此处记为 ),并构建其标准正交基 ,其中 是通过将 相对于 进行正交归一化得到的,
| (132) |
因此(参见图 6 - 右)
| (133) |
那么,我们只需在平面 上将 旋转角度 的一部分,即可得到球面插值,
| (134) |
方法三
一种类似的方法归功于 Glenn Davis 在 Shoemake(1985) 中的工作,其原理在于:连接 与 的大圆弧上任意点都必须是这两个端点的线性组合(因为这三个向量共面)。通过公式(131)计算出角度 后,我们可以从(133)中分离出 并代入(134)。应用恒等式 后,即可得到 Davis 公式(另类推导可参阅 Eberly(2010))。
| (135) |
该公式具有对称性的优势:定义反向插值器 即可得到
这与交换 和 角色后的公式完全相同。
图 7: 确保在表示 和 方向的姿态间沿最短路径进行球面线性插值。左图:四维空间中的四元数旋转平面,显示初始与最终姿态四元数,以及两种可能的插值路径:从 到 的 路径,和从 到 的 路径。右图:三维空间中的矢量旋转平面:由于 ,我们得到 ,即两个四元数产生相同的旋转。然而,插值四元数 生成的矢量 会沿从 到 的长路径旋转,而修正后的 产生 ,使得矢量 沿从 到 的最短路径旋转。
2.8 四元数与等斜旋转:揭秘其神奇之处
本节将为关于四元数的两个引人入胜的问题提供几何解释,我们称之为"魔法":
-
•
为什么乘积 会使向量 发生旋转? -
•
为何在通过 构造四元数时需要考虑半角?
首先,让我们在此重述方程(116),该方程通过四元数乘积矩阵 和 (定义见 19)表达了四元数的旋转作用。
对于单位四元数 ,其四元数乘积矩阵 和 满足两个显著特性,
| (136) | ||||
| (137) |
因此它们属于 的元素,即 空间中的正规旋转矩阵。更具体地说,它们代表了一种特殊类型的旋转,称为等斜旋转 ,我们将在下文解释。因此,根据(116),一个四元数旋转对应于 中两个串联的等斜旋转。
为了阐释四元数旋转的深层原理,我们需要理解 中的等斜旋转。为此,我们首先需要理解 中的一般旋转。而要理解 中的旋转,我们需要回溯到 ——其旋转本质上是平面旋转。让我们逐步解析这些概念。
中的旋转:
在 中,我们考虑向量 绕任意轴 的旋转——参见图 8,并回顾图 1。旋转时,平行于旋转轴 的向量保持静止,而垂直于旋转轴的向量则在垂直于该轴的平面 内旋转。对于一般向量 ,其在该平面内的两个分量会进行旋转,而轴向分量则保持静止。
图 8: 中的旋转。向量 绕轴 的旋转在垂直于该轴的平面内描绘出一个圆周。平行于轴的 分量保持不动,由轴上红色小点表示。 右侧示意图展示了平面子空间与轴子空间上旋转点的截然不同行为。
中的旋转:
在 中(参见图 9),由于额外的维度, 中的一维旋转轴变成了一个新的二维平面。
图 9: 中的旋转。可以在两个正交平面 和 上进行两种正交旋转。向量 (未绘制)在平面 中的旋转会导致该平面内两个平行分量( 上的红点)形成圆周运动(红色圆圈),而 中的另外两个分量保持不变(红点)。相反,在平面 中的旋转( 上的蓝点沿蓝色圆圈运动)会使 中的分量保持不变(蓝点)。右侧示意图通过放弃绘制 中无法表现的透视角度,更清晰地展示了这一情形,避免产生误导。
这第二个平面为第二次旋转提供了空间。
事实上,四元数空间中的旋转包含了四维空间中两个正交平面上的独立旋转。这意味着每个平面上的四维向量仅在其所属平面内旋转,且针对某一平面的旋转不会影响向量在另一正交平面上的分量。因此,这些平面被称为 “不变平面”。
中的等斜旋转:
等斜旋转(源自希腊语,iso:“相等”,klinein:“倾斜”)是指 中两个不变平面旋转角度大小相等的旋转。当两个角度符号相同时, 131313Given the two invariant planes, we arbitrarily select their orientations so that we can associate positive and negative rotation angles in them.
我们称之为左等斜旋转 ;当符号相反时,则称为右等斜旋转 。等斜旋转的一个显著特性(我们已在式 22 中见过)是:左、右等斜旋转具有可交换性。
| (138) |
与 中的四元数旋转:
给定单位四元数 (表示绕轴 旋转角度 的 空间旋转),矩阵 是 空间中对应于左乘四元数 的左等斜旋转,而 则是对应于右乘四元数 的右等斜旋转。这些等斜旋转的角度大小精确为 和 141414This can be checked by extracting the eigenvalues of the isoclinic rotation matrices: they are formed by pairs of conjugate complex numbers with a phase equal to . 。
且不变平面相同。
那么,旋转表达式(116),在此再次重现,
表示对 4-向量 进行两次链式等斜旋转,一次左旋一次右旋,每次旋转 中所需旋转角度的一半。在 的某个不变平面内(见图 10),这两个半角因符号相反而相互抵消。
图 10: 中的四元数旋转。 两个链式等斜旋转,一个左旋(具有相等的半角)和一个右旋(具有相反的半角),仅在其中一个不变平面内产生一个完整角度的纯旋转。
在另一个平面中,它们叠加产生总旋转角度 。若根据(116)定义所得旋转矩阵 ,可轻易发现(另见(138)),
| (139) |
其中 是 中的旋转矩阵,它显然旋转了 的 子空间中的向量,同时保持第四维度不变。
本文讨论的内容有些超出了当前文档的范围。
此外,这一论述并不完整,因为它除了在(139)中给出的结果外,并未提供为何需要进行 而非例如 的直观解释或几何说明。 151515Let it suffice to say that works for rotations if is a unit quaternion. In fact, the product produces reflections (not rotations!) in if is a unit pure quaternion. Finally, the product , with a unit non-pure quaternion, exhibits no remarkable properties.
我们在此处包含这一内容,仅作为提供另一种解释四元数旋转方式的手段,希望读者能对其机制有更直观的理解。
感兴趣的读者建议参阅 中关于等斜旋转的相关文献。
3 四元数约定:我的选择
3.1 四元数的不同变体
确定四元数的方法有多种,主要涉及以下四个二元选择:
-
•
其元素的排列顺序——实部在前或实部在后:(140) -
•
四元数代数的乘法公式定义:
这意味着,给定一个旋转轴𝐮,一个四元数𝐪right{𝐮θ}使用右手法则将向量绕𝐮旋转θ角度,而另一个四元数𝐪left{𝐮θ}则使用左手法则。 -
•
旋转算子的功能——旋转坐标系或旋转向量:Passive vs. Active. (142) -
•
在被动情况下,操作的方向——局部到全局或全局到局部:(143)
这种选择的多样性导致了 12 种不同的组合。历史发展使得某些约定优于其他约定 (Chou, 1992; Yazell, 2009)。如今,在现有文献中,我们可以找到多种四元数变体,如 Hamilton、STS 161616Space Transportation System, commonly known as NASA’s Space Shuttle. 、JPL 171717Jet Propulsion Laboratory. 、ISS 181818International Space Station. 、ESA 191919European Space Agency. 、工程学、机器人学等,可能还有更多命名方式。其中许多形式可能是相同的,有些则不同,但这一事实很少被明确说明,许多著作对上述四种选择中其四元数的描述也往往不够充分。
这些差异以非显而易见的方式影响着旋转、组合等操作的相应公式。因此,这些公式并不兼容,我们需要从一开始就做出明确选择。
两种最常用且文献记载最完备的约定是哈密顿约定((140–143)左侧选项)和 JPL 约定(右侧选项,(142 除外)。表 2 总结了它们的特性。JPL 约定主要应用于航空航天领域,而哈密顿约定更常见于机器人学等其他工程领域——但这并非绝对规则。
表 2: 关于 4 个二元选择的哈密顿与 JPL 四元数约定对比
| 四元数类型 | 哈密顿 | JPL | |
|---|---|---|---|
| 1 | 分量顺序 | ||
| 2 | 代数 | ||
| 手性 | 右手系 | 左手系 | |
| 3 | 函数 | 被动式 | 被动式 |
| 4 | 右乘运算表示 |
局部到全局 | 全局到局部 |
| 默认表示法, | |||
| 默认操作 |
如公式(2)所示,我选择了哈密顿约定。该约定采用右手系,与机器人学中广泛使用的诸多软件库(如 Eigen、ROS、Google Ceres)兼容,也与大量基于 IMU 的姿态估计卡尔曼滤波文献一致 (Chou,, 1992; Kuipers,, 1999; Piniés 等,, 2007; Roussillon 等,, 2011; Martinelli,, 2012,以及其他众多文献)。
JPL 约定可能较少被采用,至少在机器人领域如此。该约定在 (Trawny 和 Roumeliotis,2005) 中有详尽阐述,这部参考著作的目标与范围与本文极为相近,但完全专注于 JPL 约定。JPL 四元数被用于 JPL 文献(显然如此)以及 Li、Mourikis、Roumeliotis 及其同事的关键论文中(参见例如 (Li 和 Mourikis,2012;Li 等人,2014)),这些论文都借鉴了 Trawny 和 Roumeliotis 的文档。这些著作是我们处理视觉-惯性里程计与 SLAM 问题时的主要灵感来源——这也正是我们的研究方向。
在本节剩余部分,我们将更深入地分析这两种四元数约定。
3.1.1 四元数分量的顺序
虽然并非最根本的区别,但汉密尔顿四元数与 JPL 四元数之间最显著的差异在于其分量顺序——标量部分分别位于首位(汉密尔顿约定)或末位(JPL 约定)。这种改变的影响相当明显,不应构成理解上的重大挑战。事实上,只要保持其他三个特征不变,即使某些将四元数实部置于末尾的文献(例如 C++库 Eigen)仍被视为采用汉密尔顿约定。
为提升表述清晰度,我们采用下标 表示四元数分量,而非常用的 标记法。当调整分量顺序时, 始终代表实部,而 是否同样适用则并不明确——某些情况下可能遇到诸如 的表示(其中 为实部且位于末尾),但在通用的 表达式中,位于末尾的实部应为 。 202020See also footnote 1.
在不同约定间转换时,必须谨慎处理涉及完整 或 四元数相关矩阵的公式,因为它们的行和/或列需要交换位置。这操作本身并不困难,但可能难以察觉,因此容易出错。
关于分量顺序的两个有趣现象是:
-
•
以实部优先表示时,四元数自然可视为扩展复数,呈现为人熟知的实数+虚数形式。这种表达方式之所以令我们感到亲切,或许正源于此。 -
•
以实部置后表示时,采用向量形式表达的四元数 , 其格式与投影三维空间中的齐次向量完全等价, ,在这两种情况下, 均明确对应于三个笛卡尔坐标轴。在处理三维几何问题时,这种对应关系使得四元数与齐次向量的代数运算更加统一,特别是(但不仅限于)当齐次向量被约束在单位球面 上时。
3.1.2 四元数代数的规范定义
哈密尔顿约定定义了 ,因此,
| (144) |
而 JPL 约定则定义了 ,因此其四元数代数变为,
| (145) |
3.1.3 旋转算子的函数
我们已经了解了如何在三维空间中旋转向量。这在 (Shuster,, 1993) 中被称为主动解释,因为操作符(这会影响所有旋转操作符)主动地旋转向量,
| (147) |
另一种理解 和 对向量 作用的方式是:假设向量保持静止,而是我们观察者的视角发生了由 或 指定的旋转。这在此处被称为坐标系变换 ,并在 (Shuster,, 1993) 中称为被动解释,因为向量本身并未移动,
| (148) |
其中 和 是两个笛卡尔参考坐标系, 和 是同一向量 在这些坐标系中的表达。具体解释和规范表示法请参阅下文。
主动解释与被动解释由互为逆运算的算子所支配,即
汉密尔顿和 JPL 均采用被动约定。
方向余弦矩阵
少数作者将被动算子理解为并非旋转算子,而是一种称为方向余弦矩阵的姿态描述方式,
| (149) |
其中每个分量 表示源坐标系中 轴与目标坐标系中 轴夹角的余弦值。我们有以下恒等式,
| (150) |
3.1.4 旋转算子的方向
在被动情况下,第二个解释来源与旋转矩阵和四元数的操作方向有关,即从局部坐标系转换到全局坐标系,还是从全局转换到局部。
给定两个笛卡尔坐标系 和 ,我们将 和 分别定义为全局坐标系和局部坐标系。"全局"与"局部"是相对概念, 即 相对于 为全局坐标系,而 相对于 为局部坐标系——换言之, 是在参考坐标系 中定义的坐标系。 212121Other common denominations for the {global, local} frames are {parent, child} and {world, body}. The first one is convenient when more than two frames are involved in a system (e.g. the frames of each moving link in a humanoid robot); the second one is convenient for a solid vehicle body (e.g. a plane, a car) moving in a unique reference frame identified as the world.
我们定义 和 分别为将向量从坐标系 转换到坐标系 的四元数和旋转矩阵,其转换关系表现为:坐标系 中的向量 在坐标系 中通过四元数乘积和矩阵乘积表示为
| (151) |
相反的转换(从 到 )则通过
| (152) |
其中
| (153) |
汉密尔顿采用局部到全局作为默认规范,表示在帧 中表达的帧 ,
| (154) |
而 JPL 则采用相反的全局到局部转换,
| (155) |
请注意
| (156) |
这并不特别有用,但说明了在混合约定时很容易混淆。
还需注意,我们可以得出 的结论,但这远非一个美妙的结果,反而是极大混淆的根源,因为等式仅存在于四元数值中,而这两个四元数在公式中使用时,其含义和代表的事物是不同的。
4 扰动、导数与积分
4.1 中的加性与减性算子
图 11:S3 流形是 中的单位球面,此处以单位圆(蓝色)表示,所有单位四元数均位于其上。该流形的切空间为超平面 ,此处以直线(红色)表示。 左图 : 和 运算符将 的元素映射到 的元素,或反之。 右图 : 和 运算符将流形上的元素与切空间中的元素相关联。(同样,这些图示也适用于 流形。)
在向量空间 中,加减运算通过常规的' '加法和' '减法实现。而在 中无法直接进行此类运算,但可通过定义等效运算符来建立完整的微积分体系。
因此,我们定义在切空间元素 与 处的元素 之间的加减运算符 如下。
加运算符
“加”运算符 生成 中的一个元素 ,该元素是将 中的参考元素 与一个(通常较小的)旋转组合的结果。此旋转由向量空间 中的一个向量指定,该向量空间在参考元素 处与 流形相切,即
| (157) |
请注意,该运算符可针对 的任何表示形式进行定义。特别地,对于四元数与旋转矩阵,我们有
| (158) | ||||
| (159) |
减运算符
“减”运算符 是上述运算的逆运算。它返回 中两个元素之间的向量角度差 。该差值在参考元素 的切向量空间中表示,
| (160) |
对于四元数和旋转矩阵,其表达式为:
| (161) | ||||
| (162) |
在这两种情况下,请注意尽管向量差值 通常假设为微小量,但上述定义适用于任意大小的 (直至 流形的首次覆盖范围,即角度 以内)。
4.2 四种可能的导数定义
4.2.1 向量空间到向量空间的函数
标量与向量情形遵循经典的导数定义:给定函数 ,我们使用 将导数定义为
| (163) |
欧拉积分生成形式如下的线性表达式
4.2.2 从 到 的函数
给定函数 ,其中 ,以及一个局部小角度变化 ,我们使用 将导数定义为
| (164) | |||||
| (165) | |||||
欧拉积分生成如下形式的表达式,
4.2.3 向量空间到 的函数
对于函数 的情况,我们使用‘+’表示向量扰动,而‘ ’表示 差异。
| (166) | |||||
| (167) | |||||
欧拉积分生成如下形式的表达式,
4.2.4 从 到向量空间的函数
对于函数 的情况,我们使用' '表示 扰动,用' '表示向量差值,
| (168) | |||||
| (169) | |||||
欧拉积分生成如下形式的表达式,
4.3 旋转的有用及极其有用的雅可比矩阵
考虑将向量 绕单位轴 旋转 弧度。我们以三种等价形式表达旋转规范,即 、 和 。我们关注旋转结果相对于不同变量的雅可比矩阵。
4.3.1 关于向量的雅可比矩阵
向量 旋转后对其自身的导数计算较为简单,
| (170) |
4.3.2 关于四元数的雅可比矩阵
4.3.3 的右雅可比矩阵
让我们考虑(见图 12)一个元素 和一个旋转向量 ,使得 。当 被改变一个量 时,元素 会发生变化。在 处的 切空间中用旋转向量 表示 的变化,我们有(请参见图示,这里并非虚构)
| (175) |
也可以写作,
| (176) |
甚至
| (177) |
图 12: 右雅可比矩阵 将参数 附近的变分 映射到流形在点 处切向量空间上的变分 。
该雅可比矩阵被称为 的右雅可比矩阵,其定义为,
| (179) |
其表达式独立于所使用的参数化方法,尽管确实可以针对每种参数化形式特别表达。利用(166)式可得,
| (180) | |||||
| if using | (181) | ||||
| (182) | |||||
右雅可比矩阵及其逆矩阵可通过闭合形式计算 (Chirikjian,, 2012, 第 40 页),
| (183) | ||||
| (184) |
的右雅可比矩阵具有以下性质:对于任意 和微小 ,
| (185) | ||||
| (186) | ||||
| (187) |
4.3.4 关于旋转向量的雅可比矩阵
4.3.5 旋转组合的雅可比矩阵
4.4 扰动、不确定性与噪声
4.4.1 局部扰动
受扰动的方向 可以表示为未受扰动方向 与一个微小局部扰动 的组合。根据哈密尔顿约定,该局部扰动出现在组合积的右侧 ——为便于比较,我们同时给出其矩阵等价形式,
| (190) |
这些局部扰动 (或 )可通过指数映射从其定义在切空间中的等效向量形式 轻松获得。由此可得
| (191) |
导出局部扰动的表达式
| (192) |
4.4.2 全局扰动
考虑全局定义的扰动及其相关导数不仅是可能的,而且确实颇具意义。全局扰动出现在组合积的左侧 ,即:
| (194) |
从而得到全局扰动的表达式
| (195) |
同样,这些扰动可以在向量空间 中指定,该空间与 流形在原点处相切。
4.5 时间导数
将局部扰动表达在向量空间中,我们便能轻松推导出时间导数的表达式。
只需将 视为原始状态, 视为扰动状态,然后应用导数的定义即可。
| (196) |
对于上述情况,有
| (197) |
其中, 作为局部角扰动,对应于由 定义的局部坐标系中的角速率向量。
这些表达式显然与旋转群 框架下推导的(99)和(67)式相同。但值得注意的是,此处我们能够明确将角速率 关联到特定参考系——在本例中即由姿态 或 定义的局部坐标系。之所以实现这一点,是因为我们赋予了算子 和 精确的几何意义。从这个视角看,(200)式描述了一个参考系姿态的演化过程,其中角速率是在该参考系局部坐标系中表达的。
与全局扰动相关的时间导数可通过类似于(198)的推导得出,其结果为
| (201) |
其中
| (202) |
该式表示在全局坐标系中表达的角速率向量。方程(201)描述了当角速率在全局参考系中表达时,参考系姿态的演化过程。
4.5.1 全局到局部关系
从前文可以注意到局部角速率与全局角速率之间存在如下关系,
| (203) |
然后,后乘以共轭四元数我们得到
| (204) |
同理,考虑到对于微小 有 ,我们可得
| (205) |
这意味着,我们可以通过四元数或旋转矩阵对角速率向量 和微小角度扰动 进行坐标系转换,就像处理常规向量一样。通过设定 或 ,并注意到旋转轴向量 遵循常规变换规律,亦可得出相同结论。
| (206) |
4.5.2 四元数乘积的时间导数
我们采用常规的乘积导数公式,
| (207) |
但需注意,由于四元数乘法不满足交换律,我们必须严格遵循操作数的顺序。
这意味着 ,与标量情况不同,而是
| (208) |
4.5.3 涉及导数的其他实用表达式
可以推导出局部旋转角速度的表达式
| (209) |
以及全局旋转速率,
| (210) |
4.6 旋转速率的时间积分
四元数形式下的旋转随时间累积是通过积分与旋转速率定义相对应的微分方程实现的:局部旋转速率定义对应方程(200),全局定义对应(201)。在我们研究的场景中,角速率由局部传感器测量,因此在离散时间点 提供局部测量值 。本文仅聚焦于此情况,对应微分方程(200)为:
| (211) |
我们开发了零阶和一阶积分方法(图 13 及 14),均基于 在时间 处的泰勒展开。记 与 , 同理。该泰勒级数展开式为:
| (212) |
通过重复应用四元数导数的表达式(211)并结合 ,可以轻松求得上述 的连续导数。我们得到
| (213a) | ||||
| (213b) | ||||
| (213c) | ||||
| (213d) | ||||
为简化符号表示,我们省略了 标记,即所有乘积及 的幂次均应按四元数乘积规则进行解释。
图 13: 角速度积分近似方法:红色:真实速度。蓝色:零阶近似(自下而上:前向、中向和后向)。绿色:一阶近似。
图 14: 两个连续时间步长(灰色和黑色箭头组)的积分方案,相同时间戳的变量已按列排列。左:前向积分。中:中向和一阶积分。右:后向积分。
4.6.1 零阶积分
前向积分
后向积分
我们也可以认为,在时间段 内的恒定速度对应于该时间段结束时测得的 速度。这可以通过对 在 附近进行泰勒展开来类似地推导得出。
| (216) |
在此需要指出的是,这是处理实时到达的运动测量数据时典型的积分方法,因为积分范围对应于最后一次测量(本例中为 ,参见图 14)。为了更清晰地表达这一点,我们可以重新标记时间索引,用 代替 ,并表示为:
| (217) |
中值积分
类似地,如果将速度视为时间段 内的中值速率(这并不一定是该时间段中点处的速度),
| (218) |
我们得到,
| (219) |
4.6.2 一阶积分
此时角速度 随时间呈线性变化。其一阶导数为常数,所有更高阶导数均为零,
| (220) | ||||
| (221) |
我们可以用 和 表示中值速率 ,
| (222) |
并推导出四元数导数(213)中出现的 幂次表达式,改用更便捷的 和 表示。
| (223a) | ||||
| (223b) | ||||
| (223c) | ||||
| (223d) | ||||
将其代入四元数导数中,并在泰勒级数(212)中进行适当重排后可得:
| (224a) | ||||
| (224b) | ||||
| (224c) | ||||
| (224d) | ||||
其中在(224a)中可识别出指数级数 ,(224b)项消失,而(224d)代表我们将忽略的高阶多重项。经简化后得到(现恢复常规 表示法):
| (225) |
| (226) |
该结果与 (Trawny 和 Roumeliotis,2005) 的结论等效,但采用哈密顿约定,并使用四元数乘积形式而非矩阵乘积形式。最后,由于 (参见(33)),我们得到替代形式:
| (227) |
在此表达式中,求和的第一项是中间零阶积分器(219)。第二项是一个二阶修正项,当 与 共线时(即当旋转轴从 到 未发生变化时),该项消失。
固定旋转轴情形
让我们记 ,并称 为旋转轴。在旋转轴 恒定的情况下,我们有 ,因此,
| (228) |
这一结果实际上对于不限于 一阶导数的情况也颇具意义。事实上,若旋转轴保持不变,旋转对四元数的无穷小贡献是可交换的,即,
因此我们得到恒等式,
| (229a) | ||||
| (229b) | ||||
| (229c) | ||||
其中 表示在时间间隔 内的总旋转角度。
旋转轴变化的情况
显然,(227)式中求和项的第二部分通过 捕捉了旋转轴变化对积分方向的影响。在实际应用中,我们注意到,考虑到典型 IMU 采样时间 ,以及由于惯性作用导致的 与 近乎共线特性,这个二阶项的量级约为 或更小。更高阶的 乘积项则更为微小,因此已被忽略。
5 基于 IMU 系统的误差状态运动学
5.1 动机
我们希望建立惯性系统的误差状态运动学方程,该系统通过整合加速度计和陀螺仪的读数(含偏置和噪声)进行工作,并采用哈密顿四元数表示空间中的姿态(attitude)。
加速度计和陀螺仪读数通常来自惯性测量单元(IMU)。对 IMU 读数进行积分会形成航位推算定位系统,这类系统会随时间产生漂移。要消除漂移,需要将此类信息与 GPS 或视觉等绝对位置读数进行融合。
误差状态卡尔曼滤波器(ESKF)是我们可用于此目的的工具之一。在卡尔曼滤波范式中,ESKF 最显著的优点如下 (Madyastha 等, 2011):
-
•
方向误差状态是最小化的( 即 ,其参数数量与自由度相同),避免了与过参数化(或冗余)相关的问题,以及由此导致的协方差矩阵奇异风险,这种风险通常源于强制约束。 -
•
误差状态系统始终在原点附近运行,因此远离可能的参数奇点、万向节锁定等问题,从而确保线性化有效性始终成立。 -
•
误差状态始终很小,这意味着所有二阶乘积均可忽略不计。这使得雅可比矩阵的计算非常简便快捷。某些雅可比矩阵甚至可能是常数或等于现有状态量值。 -
•
由于所有大信号动态都已整合至标称状态中,误差动态变化缓慢。这意味着我们可以以低于预测频率的速率应用卡尔曼滤波校正(这是观测误差的唯一手段)。
5.2 误差状态卡尔曼滤波器详解
在误差状态滤波器的表述中,我们讨论真实状态、标称状态和误差状态值,其中真实状态被表示为标称状态与误差状态的适当组合(线性求和、四元数乘积或矩阵乘积)。其核心思想是将标称状态视为大信号(可进行非线性积分),而将误差状态视为小信号(因此可线性积分且适用于线性高斯滤波)。
误差状态滤波器可解释如下。
一方面,高频 IMU 数据 被积分至名义状态 。该名义状态未考虑噪声项 及其他可能的模型缺陷,因此会不断累积误差。这些误差被汇集于误差状态 中,并通过误差状态卡尔曼滤波器(ESKF)进行估计——此时所有噪声与扰动均被纳入计算。误差状态由小信号量值构成,其演化函数由(时变)线性动态系统正确定义,其动态矩阵、控制矩阵及测量矩阵均根据名义状态值计算得出。在名义状态积分的同时,ESKF 对误差状态进行高斯估计预测。此时仅执行预测步骤,因为尚无其他测量数据可用于修正估计值。当非 IMU 的观测信息( 例如 GPS、视觉等)到达时——这些信息能使误差具有可观测性且通常以远低于积分阶段的频率出现——滤波器执行修正步骤,从而获得误差状态的后验高斯估计。 此后,将误差状态的均值注入名义状态,然后将其重置为零。误差状态的协方差矩阵会相应更新以反映这一重置操作。系统将如此持续运行。
5.3 连续时间系统运动学
所有相关变量的定义总结于表 3 中。关于约定有两个重要决策值得注意:
-
•
角速度 被定义为相对于名义四元数的局部量 。这使得我们可以直接使用陀螺仪测量值 ,因为它们提供的是机体参考系下的角速度。 - •
表 3: 误差状态卡尔曼滤波器中的所有变量。
| 幅值 | 真实值 | 标称值 | 误差 | 复合运算 | 测量值 | 噪声 |
|---|---|---|---|---|---|---|
完整状态(1) |
||||||
| 位置 | ||||||
| 速度 | ||||||
四元数(2,3) |
||||||
旋转矩阵(2,3) |
||||||
角度向量(4) |
||||||
| 加速度计偏置 | ||||||
| 陀螺仪偏置 | ||||||
| 重力矢量 | ||||||
| 加速度 | ||||||
| 角速率 | ||||||
(1) 符号 表示通用组合关系 |
||||||
(2) 表示非最小化表达形式 |
||||||
(3) 全局定义角度误差情况下的组合公式参见表 4 |
||||||
(4) 指数定义如(101)及(69, 77)式所示 |
||||||
5.3.1 真实状态运动学
真实运动学方程为
| (230a) | ||||
| (230b) | ||||
| (230c) | ||||
| (230d) | ||||
| (230e) | ||||
| (230f) | ||||
此处,真实加速度 和角速度 由 IMU 以带噪声的传感器读数 和 的形式在机体坐标系中获取,即 232323It is common practice to neglect the Earth’s rotation rate in the rotational kinematics described in (232), which would otherwise be .
Considering a non-null Earth rotation rate is, in the vast majority of practical cases, unjustifiably complicated.
However, we notice that when employing high-end IMU sensors with very small noises and biases, a value of /h rad/s might become directly measurable; in such cases, in order to keep the IMU error model valid, the rate should not be neglected in the formulation.
| (231) | ||||
| (232) |
通过 。由此,真实值可以被分离出来(这意味着我们已经反转了测量方程),
| (233) | ||||
| (234) |
将上述内容代入后得到运动学系统
| (235a) | ||||
| (235b) | ||||
| (235c) | ||||
| (235d) | ||||
| (235e) | ||||
| (235f) | ||||
我们可以将其命名为 。该系统具有状态 ,由 IMU 噪声读数 控制,并受到白高斯噪声 的扰动,所有这些都由以下定义:
| (236) |
需注意的是,在上述公式体系中,重力矢量@g0#将由滤波器进行估计。该矢量遵循恒定演化方程(235f),这与已知恒定的物理量特性相符。系统始于固定且任意设定的初始姿态@g1#,由于该姿态通常不在水平面内,导致初始重力矢量通常未知。为简化起见,通常取@g2#,因此@g3#。我们估计的是在@g5#坐标系中表达的重力@g4#,而非水平坐标系中的@g6#,从而将初始姿态不确定性转化为重力方向的初始不确定性。这种处理旨在提升线性度:实际上,方程(235b)现在对承载所有不确定性的@g7#呈线性关系,而初始姿态@g8#被视为确知量,因此@g9#的初始状态不存在不确定性。一旦重力矢量完成估计,即可重建水平面,若需要还可将整个状态量及重建的运动轨迹重新定向以反映估计的水平基准。 详见 (Lupton and Sukkarieh,, 2009) 的进一步论证。这当然是可选的,读者可以自由移除系统中所有与重力相关的方程,转而采用更经典的方法,即考虑 ,其中 表示实验现场重力向量的有效小数位数,以及不确定的初始方向 。
5.3.2 标称状态运动学
标称状态运动学对应于无噪声或扰动的建模系统,
| (237a) | ||||
| (237b) | ||||
| (237c) | ||||
| (237d) | ||||
| (237e) | ||||
| (237f) | ||||
5.3.3 误差状态运动学
目标是确定误差状态的线性化动态特性。针对每个状态方程,我们写出其组合形式(见表 3),求解误差状态并简化所有二阶无穷小量。下文给出完整的误差状态动态系统,随后进行注释和证明。
| (238a) | ||||
| (238b) | ||||
| (238c) | ||||
| (238d) | ||||
| (238e) | ||||
| (238f) | ||||
方程(238a)、(238d)、(238e)和(238f)分别描述了位置误差、两种偏差误差以及重力误差,这些方程均由线性方程推导而来,其误差状态动力学较为简单。例如,考虑真实位置方程(235a)与标称位置方程(237a),根据表 3 中的组合关系 ,求解 即可得到方程(238a)。
方程(238b):线性速度误差
我们需要确定速度误差的动力学关系 。首先建立以下关系式:
| (239) | ||||
| (240) |
其中(239)是 的小信号近似,而在(240)中我们重写了(237b),但引入了 和 ,定义为体坐标系下的大信号和小信号加速度,
| (241) | ||||
| (242) |
因此我们可以将惯性系中的真实加速度表示为大信号项与小信号项的叠加,
| (243) |
我们通过以两种不同形式(左侧和右侧展开)书写 的表达式(235b)来继续推导,其中忽略了 项,
在左右两边消去 后,这导致
| (244) |
消除二阶项并重新组织一些交叉乘积(使用 ),我们得到
| (245) |
| (246) |
经适当整理后可得线速度误差的动力学方程,
| (247) |
为了进一步简化该表达式,我们通常可以假设加速度计噪声是白噪声、不相关且各向同性的 242424This assumption cannot be made in cases where the three accelerometers are not identical. ,
| (248) |
即协方差椭球是一个以原点为中心的球体,这意味着其均值和协方差矩阵在旋转下保持不变( 证明: 和 )。因此我们可以重新定义加速度计噪声向量,且完全不影响结果,根据
| (249) |
从而得到
| (250) |
方程(238c):姿态误差。
我们希望确定角误差的动态特性 。我们从以下关系出发
| (251) | ||||
| (252) |
这些就是四元数导数的真实定义与标称定义。
如同处理加速度时那样,为清晰起见,我们将角速率中的大信号项与小信号项分组,
| (253) | ||||
| (254) |
因此 可以表示为标称部分与误差部分之和
| (255) |
5.4 离散时间系统运动学
上述微分方程需要被积分成差分方程以适应离散时间间隔 。积分方法可能各不相同:某些情况下可采用精确的闭式解,其他情况下则可能使用不同精度的数值积分方法。具体积分方法详见附录。
需对以下子系统进行积分:
-
1.
标称状态。
-
2.
误差状态
-
(a)
确定性部分:状态动力学与控制。 -
(b)
随机部分:噪声与扰动。
-
(a)
5.4.1 标称状态运动学
我们可以将名义状态差分方程表示为
| (260a) | ||||
| (260b) | ||||
| (260c) | ||||
| (260d) | ||||
| (260e) | ||||
| (260f) | ||||
其中 表示 类型的时间更新, 是与当前标称姿态 相关联的旋转矩阵,而根据(101), 是与旋转 对应的四元数。
我们也可以采用更精确的积分方法,更多信息请参阅附录。
5.4.2 误差状态运动学
确定性部分按常规方式积分(此处我们遵循附录 C.2 中的方法),而随机部分的积分则产生随机脉冲(参见附录 E),因此,
| (261a) | ||||
| (261b) | ||||
| (261c) | ||||
| (261d) | ||||
| (261e) | ||||
| (261f) | ||||
此处, 、 、 和 表示作用于速度、姿态及偏差估计的随机脉冲,由高斯白噪声过程建模。其均值为零,协方差矩阵通过对 、 、 和 的协方差在步长时间 内积分获得(参见附录 E)。
| (262) | |||||
| (263) | |||||
| (264) | |||||
| (265) |
其中 、 、 和 需根据 IMU 数据手册信息或实验测量确定。
5.4.3 误差状态雅可比矩阵与扰动矩阵
雅可比矩阵可通过直接观察前一节中的误差状态差分方程获得。
为以紧凑形式表示这些方程,我们定义标称状态向量 、误差状态向量 、输入向量 及扰动脉冲向量 如下(具体推导与论证参见附录 E.1):
| (266) |
此时误差状态系统为
| (267) |
ESKF 预测方程可表示为:
| (268) | ||||
| (269) |
式中 252525 means that is a Gaussian random variable with mean and covariances matrix specified by and . ; 和 表示 对误差向量及扰动向量的雅可比矩阵; 为扰动脉冲的协方差矩阵。
上述雅可比矩阵和协方差矩阵的表达式将在下文详细说明。此处出现的所有状态相关值均直接取自标称状态。
| (270) |
| (271) |
最后请注意,你绝对不应跳过协方差预测(269)!!实际上,项 非零,因此该协方差会持续增长——这正是任何预测步骤中应有的表现。
6 IMU 与互补传感器数据的融合
当接收到除 IMU 外的其他信息(如 GPS 或视觉数据)时,我们着手修正 ESKF。在精心设计的系统中,这将使 IMU 偏差可观测,并允许 ESKF 准确估计这些偏差。存在多种组合可能,最流行的包括 GPS+IMU、单目视觉+IMU 以及立体视觉+IMU。近年来,视觉传感器与 IMU 的组合引起了广泛关注,并因此催生了大量科研活动。这些视觉+IMU 配置在 GPS 拒止环境中极具应用价值,既可部署于移动设备(通常是智能手机),也可应用于无人机等小型敏捷平台。
虽然迄今为止 IMU 信息一直用于 ESKF 的预测环节,但其他类型的信息则用于修正滤波器,从而观测 IMU 的偏差误差。修正过程包含三个步骤:
-
1.
通过滤波器修正观测误差状态, -
2.
将观测误差注入标称状态,且 -
3.
重置误差状态。
这些步骤将在后续章节中展开说明。
6.1 通过滤波器校正观测误差状态
假设我们有一个传感器,它提供依赖于状态的信息,例如
| (272) |
其中 是系统状态(真实状态)的一般非线性函数, 是具有协方差 的高斯白噪声,
| (273) |
我们的滤波器正在估计误差状态,因此滤波器校正方程 262626We give the simplest form of the covariance update, .
This form is known to have poor numerical stability, as its outcome is not guaranteed to be symmetric nor positive definite.
The reader is free to use more stable forms such as 1) the symmetric form and 2) the symmetric and positive Joseph form . ,
| (274) | ||||
| (275) | ||||
| (276) |
需要定义关于误差状态 的雅可比矩阵 ,并在最佳真实状态估计 处求值。由于此阶段误差状态均值为零(我们尚未观测到它),我们有 ,因此可以使用标称误差 作为求值点,从而得到
| (277) |
6.1.1 滤波器校正的雅可比矩阵计算
上述雅可比矩阵可通过多种方式计算。最具启发性的方法是利用链式法则:
| (278) |
此处 是 关于其自身参数的标准雅可比矩阵( 即常规 EKF 中使用的雅可比矩阵)。链式法则的第一部分取决于所用特定传感器的测量函数,此处不作展示。
第二部分 是真实状态相对于误差状态的雅可比矩阵。这部分可以在此推导,因为它仅依赖于误差状态卡尔曼滤波器(ESKF)的状态组合。我们得到以下导数:
| (279) |
这将导致除四元数项 (例如 )之外的所有 块均为单位矩阵(如 所示)。因此我们得到如下形式:
| (280) |
6.2 观测误差注入标称状态
在 ESKF 更新后,名义状态通过适当的组合运算(加法或四元数乘法,参见表 3)与观测到的误差状态进行更新:
| (282) |
即,
| (283a) | ||||
| (283b) | ||||
| (283c) | ||||
| (283d) | ||||
| (283e) | ||||
| (283f) | ||||
6.3 误差状态卡尔曼滤波重置
将误差注入名义状态后,误差状态均值 将被重置。这对于姿态部分尤为重要,因为新的姿态误差将相对于新名义状态的姿态坐标系进行局部表达。为使 ESKF 更新完整,误差协方差需根据此调整进行相应更新。
我们称误差重置函数为 ,其表达式如下:
| (284) |
其中 代表 的复合逆。因此,ESKF 误差重置操作可表示为:
| (285) | ||||
| (286) |
此处 为由下式定义的雅可比矩阵:
| (287) |
与前述更新雅可比矩阵的情况类似,该雅可比矩阵在所有对角块上均为单位矩阵,仅方向误差块例外。
此处给出完整表达式,并在下一节中推导方向误差块 。
| (288) |
在大多数情况下,误差项 可被忽略,仅需计算雅可比矩阵 ,从而实现简单的误差重置。这是误差状态卡尔曼滤波器(ESKF)多数实现采用的方式。本文提供的表达式应能产生更精确的结果,对于减少里程计系统的长期误差漂移可能具有重要意义。
6.3.1 针对姿态误差的重置操作雅可比矩阵
我们需要推导新角误差 相对于旧误差 和观测误差 的表达式。考虑以下事实:
-
•
真实方位在误差重置时不会改变, 即 。由此可得:(289) - •
结合这两个恒等式,我们得到了 的表达式,
| (291) |
考虑到
上述恒等式可以展开为
| (292) |
由此得到一个标量方程和一个向量方程,
| (293a) | ||||
| (293b) | ||||
其中第一个关系式信息量不大,因为它仅表达了无穷小量之间的关系。
从第二个方程可以证明 ,这正是重置操作所预期的结果。通过简单观察即可得到雅可比矩阵,
| (294) |
7 采用全局角度误差的误差状态卡尔曼滤波器
角度误差的全局定义 意味着左侧合成 , 即 ,
出于完备性考虑需要说明的是:无论角度误差是否采用全局定义,我们仍保持角速度矢量 的局部定义, 即连续时间下的 ,以及离散时间下的 。这样处理是出于便利性考虑,因为陀螺仪提供的角速度测量值始终基于载体坐标系(即局部坐标系)。
7.1 连续时间下的系统运动学
7.1.1 真实状态与标称状态的运动学
真实运动学与标称运动学不涉及误差,其方程保持不变。
7.1.2 误差状态运动学
我们首先写出误差状态运动学的方程,随后进行相关注释与证明。
| (295a) | ||||
| (295b) | ||||
| (295c) | ||||
| (295d) | ||||
| (295e) | ||||
| (295f) | ||||
其中,除 和 外的所有方程均为平凡表达式。非平凡表达式将在下文中展开推导。
方程(295b):线速度误差
方程(295c):方向误差。
我们首先写出四元数导数的真实定义与名义定义,
| (301) | ||||
| (302) |
需特别说明的是,此处采用全局定义的角误差, 即
| (303) |
如同处理局部定义误差的情况一样,我们也将大信号与小信号角速率(253–254)进行分组。我们通过两种不同方式(左展开与右展开)计算 ,
得到 后,可简化为
| (304) |
将左右项右乘 ,并回顾 ,可进一步推导如下,
| (305) |
其中 表示全局坐标系下表达的小信号角速率。接着,
| (306) |
最终导出一个标量等式和一个向量等式
| (307a) | ||||
| (307b) | ||||
第一个方程得出 ,它由二阶无穷小量构成,实用性不强。第二个方程在忽略所有二阶项后得到,
| (308) |
最后,回顾(254),我们得到了全局角度误差的线性化动力学方程,
| (309) |
7.2 离散时间下的系统运动学
7.2.1 标称状态
标称状态方程不涉及误差项,因此与姿态误差在局部定义时的情况相同。
7.2.2 误差状态
采用欧拉积分法,我们得到以下差分方程组,
| (310a) | ||||
| (310b) | ||||
| (310c) | ||||
| (310d) | ||||
| (310e) | ||||
| (310f) | ||||
7.2.3 误差状态雅可比矩阵与扰动矩阵
通过直接观察上述方程即可获得状态转移矩阵,
| (311) |
在考虑各向同性噪声及附录 E 的推导后,扰动雅可比矩阵与扰动矩阵保持不变,
| (312) |
7.3 与互补传感器数据的融合
当考虑全局角度误差时,涉及误差状态卡尔曼滤波(ESKF)机制的融合方程仅需微调。我们通过 ESKF 校正环节修订了误差状态观测中的这些变化,包括将误差注入标称状态的过程以及重置步骤。
7.3.1 误差状态观测
与局部误差定义唯一的不同之处在于观测函数中涉及姿态与角度误差关系的雅可比矩阵块。下文将详细推导这一新的矩阵块。
7.3.2 观测误差向标称状态的注入
7.3.3 ESKF 重置
根据以下公式重置 ESKF 误差均值并更新协方差矩阵,
| (315) | ||||
| (316) |
其雅可比矩阵为
| (317) |
其非平凡项推导如下。
我们的目标是获得新角误差 相对于旧误差 的表达式。我们考虑以下事实:
-
•
真实方位在误差重置时不会改变, 即 。由此可得:(318) - •
结合这两个恒等式,我们得到了新姿态误差相对于旧误差及观测误差 的表达式。
| (320) |
考虑到
上述恒等式可以展开为
| (321) |
由此得到一个标量方程和一个向量方程,
| (322a) | ||||
| (322b) | ||||
其中第一个关系式信息量不大,因为它仅表达了无穷小量之间的关系。
从第二个方程可以证明 ,这正是重置操作所预期的结果。通过简单观察即可得到雅可比矩阵,
| (323) |
与局部误差情况的差异总结在表 4 中。
表 4: 与方向误差定义相关的算法修改。
| 上下文 | 项目 | 局部角度误差 | 全局角度误差 |
|---|---|---|---|
| 误差组合 | |||
| 欧拉积分 | |||
| 误差观测 | |||
| 误差注入 | |||
| 误差重置 |
附录 A 龙格-库塔数值积分方法
我们的目标是积分如下形式的非线性微分方程
| (324) |
在有限时间间隔 内,为了将其转换为差分方程, 即 ,
| (325) |
或者说,如果我们假设 和 ,
| (326) |
最常用的方法族之一是龙格-库塔方法(以下简称 RK)。
这些方法通过多次迭代来估计区间内的导数,然后利用该导数进行步长 的积分。
在接下来的章节中,将介绍几种龙格-库塔方法,从最简单到最通用的形式,并按其最常用的名称进行命名。
注:此处所有材料均摘自英文维基百科的龙格-库塔方法条目。
A.1 欧拉方法
欧拉方法假设导数 在区间内保持恒定,因此
| (327) |
作为一种通用的龙格-库塔方法,这对应于单阶段方法,其计算过程可描述如下:首先在初始点处计算导数,
| (328) |
并利用它计算终点处的积分值,
| (329) |
A.2 中点方法
中点法假设导数为区间中点处的导数,并通过一次迭代计算 在该中点处的值, 即
| (330) |
中点法可以解释为如下两步法。
首先,使用欧拉法积分至中点,采用先前定义的 ,
| (331) | ||||
| (332) |
然后利用该值计算中点处的导数 ,从而完成积分
| (333) | ||||
| (334) |
A.3 龙格-库塔四阶方法
该方法通常简称为龙格-库塔法。
该方法假设在区间的起点、中点和终点处对 进行评估取值,并通过四个阶段或迭代计算积分,依次获得四个导数 。这些导数或斜率随后进行加权平均,得到区间内导数的四阶估计值。
相较于前述两种单步公式方法,RK4 方法更适合被描述为一个小型算法。RK4 积分步骤如下:
| (335) |
即增量是通过假设一个斜率为各斜率 的加权平均来计算,
| (336) | ||||
| (337) | ||||
| (338) | ||||
| (339) |
不同斜率的解释如下:
-
•
为区间起点处的斜率,采用 (欧拉法); -
•
为区间中点处的斜率,采用 (中点法); -
•
再次表示中点处的斜率,但此时使用的是 ; -
•
表示区间终点处的斜率,使用的是 。
A.4 通用龙格-库塔方法
更精细的龙格-库塔方法是可能实现的,其目标在于减小误差和/或增强稳定性。这些方法的一般形式为
| (340) |
其中
| (341) |
即迭代次数(方法的阶数)为 ,平均权重由 定义,计算时间点由 确定,而斜率 则通过值 计算得出。根据项 的结构不同,可以得到显式或隐式的 RK 方法。
-
•
在显式方法中,所有 都是按顺序计算的, 即仅使用先前计算的值。这意味着矩阵 是一个对角元为零的下三角矩阵( 即 对于 )。欧拉法、中点法和 RK4 法都属于显式方法。 -
•
隐式方法具有完整的 矩阵,需要通过求解线性方程组来确定所有 。因此其计算成本更高,但相比显式方法能够提高精度和稳定性。
更多详细信息请参阅专业文档。
附录 B 闭式积分方法
在许多情况下,积分步骤可以得到闭式表达式。我们现在考虑一阶线性微分方程的情况,
| (342) |
也就是说,该关系在区间内是线性且恒定的。
在此类情况下,在区间 上的积分将得到
| (343) |
其中, 被称为转移矩阵。该转移矩阵的泰勒展开式为
| (344) |
在针对已知的 实例展开此级数时,有时可以识别结果中的已知级数。这使得能够以闭合形式写出最终的积分结果。以下列举几个示例。
B.1 角度误差的积分
B.2 简化 IMU 示例
考虑由误差状态动力学支配的简化 IMU 驱动系统,
| (350a) | ||||
| (350b) | ||||
| (350c) | ||||
其中 代表 IMU 读数,我们忽略了重力和传感器偏差。该系统由状态向量和动态矩阵定义,
| (351) |
随着
| (352) | ||||
| (353) | ||||
| (354) |
我们可以观察到, 的递增幂次项包含固定部分和 的递增幂次。这些幂次可能导出封闭形式的解,如前一节所示。
我们将矩阵 划分如下,
| (357) |
让我们逐步推进,逐一探索 的所有非零分块。
平凡对角项
从对角线上的前两项开始,如所示,它们是单位矩阵。
旋转对角项
接下来是旋转对角项 ,它将新的角度误差与旧的角度误差相关联。对该项进行完整的泰勒级数展开可得
| (358) |
正如我们在前一节所见,这对应于我们熟知的旋转矩阵,
| (359) |
位置-速度项
最简单的非对角项是 ,其表达式为
| (360) |
速度-角度项
现在让我们转向项 ,通过写出其级数展开式,
| (361) |
| (362) |
这简化为
| (363) |
此时我们有两种选择。
我们可以将级数截断至第一显著项,得到 ,但这将不是一个闭合形式。下一节将展示使用这种简化方法的结果。或者,我们可以提取公因子 并重写为
| (364) |
随着
| (365) |
该级数 与我们为 所写的级数(358)相似,但存在两点差异:
-
•
在 中的幂次与有理系数 不匹配,且与 的幂次不一致。 事实上,我们在此指出, 中的下标“1”表示每个成员中都缺少一个 的幂次。 -
•
序列开头缺少了一些项。同样,下标“1”表明缺少了这样一项。
位置-角度交叉项
让我们最终处理 项。其泰勒展开式为,
| (370) |
我们提取出常数项后得到,
| (371) |
随着
| (372) |
此处我们注意到 中的下标"2"具有如下解释:
-
•
该级数的每一项都缺少 的两次幂 -
•
该级数的前两项缺失。
B.3 完整 IMU 示例
为了将简化 IMU 示例中展示的方法推广到一般情况,我们需要更细致地考察完整 IMU 案例。
| (378) |
如前所述,让我们写出 的若干幂次,
基本上,我们观察到以下现象,
-
•
对角线上唯一的旋转项 ,在幂次序列 中向右上方传播。所有未受此传播影响的项均消失。这种传播在以下三个方面影响了序列 的结构: -
•
的幂次稀疏性在第3次幂后趋于稳定。也就是说,对于高于3次的 幂次,不再有新的非零块出现或消失。 -
•
左上角的 块对应于前例中的简化 IMU 模型,与该示例相比未发生变化。因此,之前推导的闭式解仍然适用。 -
•
与陀螺仪偏差误差相关的项(第五列的那些项)引入了类似的 幂级数,可采用简化示例中使用的相同技术求解。
此刻我们关注的是寻找一种通用方法,来构建转移矩阵 闭式元素的框架。让我们回顾先前关于级数 和 的论述。
-
•
下标对应于级数各项中 的缺失幂次。 -
•
下标与级数起始处缺失的项数一致。
基于这些性质,我们引入由 272727Note that, being a square matrix that is not necessarily invertible (as it is the case for ), we are not allowed to rearrange the definition of with . 定义的级数 ,
| (379) |
其中求和从 项开始,且各项缺少矩阵 的 次幂。由此立即得出 和 满足
| (380) |
且 。现在可将转移矩阵(377)表示为这些级数的函数,
| (381) |
我们的问题现已转化为寻找 的通用闭式表达式。让我们观察目前已获得的闭式结果,
| (382) | ||||
| (383) | ||||
| (384) |
为了推导 ,我们需要两次应用恒等式(366)(因为我们缺少三次幂,而每次应用(366)仅能使该数值增加二),从而得到
| (385) |
这将导出
| (386) |
通过仔细检查级数 ,我们现在可以推导出 的通用闭式表达式,如下所示,
| (387a) | |||||
| (387b) | |||||
| (387c) | |||||
| (387d) | |||||
将 的适当值代入(381)的对应位置后,转移矩阵 的最终结果便立即得出。
值得注意的是,这些出现在 新表达式中的级数现在具有有限项数,因此可以实际计算。也就是说,只要 (这总是成立), 的表达式就是闭合形式。在当前示例中,如(381)所示,我们有 。
附录 C 基于截断级数的近似计算方法
在前文章节中,我们已推导出以线性化误差状态形式表达的复杂 IMU 驱动动态系统状态转移矩阵的闭式解。虽然闭式解始终具有研究价值,但高阶误差对实际算法性能的影响程度尚不明确。这一观点对于 IMU 积分误差以较高频率被观测(并补偿)的系统尤为重要,例如视觉-惯性或 GPS-惯性融合方案。
本节我们设计了近似求解转移矩阵的方法。这些方法基于相同的假设前提——转移矩阵可表示为泰勒级数展开式,随后在最具显著性的项处截断该级数。截断操作既可针对整个系统实施,亦可分块进行。
C.1 系统级截断
C.1.1 一阶截断:有限差分法
一种典型且广泛适用于此类系统的积分方法是
该方法基于有限差分法计算导数,即
| (388) |
这直接导致了
| (389) |
这正是欧拉积分方法的表述。
在积分区间开始时对函数 进行线性化处理会导致
| (390a) | |||
| where is a Jacobian matrix. This is strictly equivalent to writing the exponential solution to the linearized differential equation and truncating the series at the linear term (i.e., the following relation is identical to the previous one), | |||
| (390b) | |||
这意味着欧拉方法(附录 A.1)、有限差分法和一阶系统泰勒截断方法都是相同的。我们得到近似转移矩阵,
| (391) |
C.1.2 N 阶截断
采用更高阶截断将提升近似转移矩阵的精度。特别值得关注的是能最大限度利用结果稀疏性的截断阶数,即超过该阶数后不再出现新非零项的情况。
对于第 B.2 节的简化 IMU 示例,该阶数为 2,结果为
| (394) |
C.2 分块截断
将转移矩阵各分块的泰勒级数截断至首项显著项,即可获得与前述闭合形式相当接近的近似解。
这意味着,不同于前文对 整体幂级数的截断处理,我们需对各分块单独考量。因此,截断分析需基于分块逐一进行。以下通过两个示例加以说明。
对于第 B.2 节的简化 IMU 示例,我们得到了级数 和 ,可将其截断如下
| (396) | |||||||||
| (397) | |||||||||
由此得到近似转移矩阵
| (398) |
该结果比前述系统全局一阶截断更精确(因为新增了右上角项),但仍易于推导计算,特别是与第 B 节给出的闭合形式相比。需注意我们已对最低阶项( 即 )采用了闭合形式。
一般情况下,除 外,其余 均可仅取其级数首项近似, 即
| (399) |
附录 D 基于龙格-库塔积分的状态转移矩阵
另一种近似求解状态转移矩阵的方法是采用龙格-库塔积分法。
当动态矩阵 在积分区间内不能视为常数时,这种方法可能是必要的, 即 ,
| (403) |
让我们重写连续时间和离散时间下定义同一系统的两个关系式。
它们涉及动态矩阵 和转移矩阵 ,
| (404) | ||||
| (405) |
这些方程允许我们通过以下两种方式展开 (左展开与右展开,请注意小圆点表示时间导数),
| (406) |
此处(406)源于注意到 是采样值。因此,
| (407) |
这与方程(404)是相同的常微分方程,只是现在应用于转移矩阵而非状态向量。需注意,由于恒等式 的存在,区间起始处的转移矩阵 始终为单位矩阵。
| (408) |
使用 RK4 与 时,我们有
| (409) |
随着
| (410) | ||||
| (411) | ||||
| (412) | ||||
| (413) |
D.1 误差状态示例
让我们考虑针对非线性时变系统的误差状态卡尔曼滤波器
| (414) |
其中 表示真实状态, 为控制输入。该真实状态是由标称状态 与误差状态 通过 运算构成的复合状态,
| (415) |
其误差状态动力学呈现线性形式,且随时间变化取决于标称状态 与控制量 , 即 ,
| (416) |
也就是说,式(403)中的误差状态动态矩阵具有 的形式。误差状态转移矩阵的动力学可表述为,
| (417) |
为了对此方程进行龙格-库塔积分,我们需要在龙格-库塔评估点处的 和 值,对于 RK4 方法这些点是 。从简单的部分开始,评估点处的控制输入 可通过当前与上次测量值的线性插值获得,
| (418) | ||||
| (419) | ||||
| (420) |
名义状态动力学应采用最佳可行的积分方法进行计算。例如,使用四阶龙格-库塔(RK4)积分法时,
由此得到评估点处的估计值,
| (421) | ||||
| (422) | ||||
| (423) |
我们注意到此处 采用了与控制信号相同的线性插值方法。考虑到龙格-库塔更新的线性特性,这并不令人意外。
无论我们通过何种方式获得标称状态值,现在都可以计算用于转移矩阵积分的 RK4 矩阵,
最终得到,
| (424) |
附录 E 随机噪声与扰动的积分方法
我们现在的目标是提供适当的方法来对动态系统中的随机变量进行积分。当然,我们无法对未知的随机值进行积分,但为了不确定性传播,我们可以对其方差和协方差进行积分。这对于建立连续特性系统(在连续时间中指定)但以离散方式估计的估计器中的协方差矩阵是必要的。
考虑连续时间动态系统,
| (425) |
其中 为状态向量, 是包含噪声 的控制信号向量,因此控制测量值为 ,而 为随机扰动向量。噪声与扰动均假定为白高斯过程,其特性由以下参数确定:
| (426) |
其中上标 表示连续时间不确定性规范,这是我们希望进行积分的对象。
控制信号中的噪声水平 与随机扰动 在本质上存在重要差异:
-
•
在离散化过程中,控制信号在时间点 处采样,具有 。测量部分显然被视为在积分区间内保持恒定, 即 ,因此采样时刻 的噪声水平也保持恒定,(427) -
•
扰动 从不被采样。
因此,这两个随机过程在 上的积分结果存在差异。让我们具体分析这一现象。
连续时间误差状态动态方程(425)可线性化为
| (428) |
随着
| (429) |
并在采样周期 内进行积分,得到
| (430) | ||||
| (431) |
该表达式包含三个性质迥异的项。可按如下方式进行积分:
-
1.
由附录 B 可知,动态部分经积分得到状态转移矩阵
其中Φ=e𝐀ΔtΦ=e^{{\bf A}\Delta t}可通过闭式解或不同精度等级的近似方法计算得出。 - 2.
-
3.
根据概率论可知,连续高斯白噪声在时间间隔Δt 内的积分将产生离散高斯白噪声脉冲𝐰𝑛。我们注意到,与前述测量噪声不同,该扰动在积分区间内不具有确定性行为,因此必须采用随机积分方法处理。
因此,离散时间误差状态动态系统可表述为
| (435) |
其状态转移矩阵、控制矩阵及扰动矩阵分别由下式给出
| (436) |
其中噪声与扰动水平定义为
| (437) |
随着
| (438) |
表 5: 积分对系统矩阵及协方差矩阵的影响。
| 描述 | 连续时间 | 离散时间 |
|---|---|---|
| 状态 | ||
| 误差状态 | ||
| 系统矩阵 | ||
| 控制矩阵 | ||
| 扰动矩阵 | ||
| 控制协方差 | ||
| 扰动协方差 |
这些结果总结在表 5 中。扩展卡尔曼滤波的预测阶段将根据以下公式传播误差状态的均值与协方差矩阵:
| (439) | ||||
| (440) |
此处观察积分区间 对协方差更新(440)三项的不同影响具有重要启示意义:动态误差项呈指数变化,测量误差项呈二次变化,而扰动误差项呈线性变化。
E.1 噪声与扰动脉冲
人们经常会遇到(例如在重用现有代码或解释其他作者的文档时)比我们这里使用的形式更简单的 EKF 预测方程,即:
| (441) |
这对应于一般的离散时间动态系统,
| (442) |
其中
| (443) |
是一个由随机(白噪声、高斯分布)脉冲组成的向量,这些脉冲在时刻 直接叠加到状态向量上。矩阵 被视为脉冲的协方差矩阵。根据我们之前的分析,应按下式计算该协方差矩阵:
| (444) |
当脉冲不作用于完整状态时(这种情况很常见),矩阵 并非全对角矩阵且可能包含大量零元素。此时可写出等效形式:
| (445) |
随着
| (446) |
其中矩阵 仅将每个独立脉冲映射至其影响的状态向量部分。相应的协方差矩阵 则更小且为完全对角矩阵。具体示例请参阅下一节。在此情况下,误差状态卡尔曼滤波器(ESKF)的时间更新变为
| (447) | ||||
| (448) |
显然,所有这些形式都是等价的,这可以从以下关于一般扰动 的双重恒等式中得到验证
| (449) |
E.2 完整 IMU 示例
我们研究了一种用于惯性测量单元(IMU)的误差状态卡尔曼滤波器的构建方法。
误差状态系统定义在(238)中,包含名义状态 、误差状态 、含噪声的控制信号 以及由以下公式指定的扰动 :
| (450) |
在本文所考虑的 IMU 模型中,控制噪声对应于 IMU 测量中的加性噪声。这些扰动会影响偏置,从而产生其随机游走行为。动态矩阵、控制矩阵和扰动矩阵分别为(参见(428)、(378)和(238)):
| (451) |
对于常规配置的三轴同类型加速度计与陀螺仪组合的 IMU 而言,噪声与扰动具有各向同性特性。其标准差按以下标量形式给出:
| (452) |
其协方差矩阵为纯对角矩阵,即
| (453) |
E.2.1 噪声与扰动脉冲
在采用脉冲形式的扰动规范 的情况下,我们可以如下重新定义系统:
| (454) |
其中名义状态、误差状态、控制及冲量向量分别定义为:
| (455) |
其转移矩阵和扰动矩阵定义为:
| (456) |
以及由以下方差指定的脉冲
| (457) |
考虑到(451)中 的定义, 的平凡规范可能显得尤为令人惊讶。实际情况是,误差在 中被定义为各向同性的,因此 和 成立,从而导出了 的表达式。但在考虑非各向同性 IMU 时,这种做法不可行,此时需要配合使用恰当的雅可比矩阵 与正确的 规范。
当然,我们也可以使用全状态扰动脉冲,
| (458) |
随着
| (459) |
再见
References
- Chirikjian, (2012) Chirikjian, G. S. (2012). Stochastic Models, Information Theory, and Lie Groups, volume 2: Analytic Methods and Modern Applications of Applied and Numerical Harmonic Analysis. Birkhäuser,, Basel.
- Chou, (1992) Chou, J. (1992). Quaternion kinematic and dynamic differential equations. Robotics and Automation, IEEE Transactions on, 8(1):53–64.
- Dam et al., (1998) Dam, E. B., Koch, M., and Lillholm, M. (1998). Quaternions, interpolation and animation. Technical report, Department of Computer Science, University of Copenhagen, Denmark.
- Eberly, (2010) Eberly, D. (2010). Quaternion algebra and calculus. Technical report, Geometric Tools, LLC.
- Kuipers, (1999) Kuipers, J. B. (1999). Quaternions and rotation sequences: a primer with applications to orbits, aerospace, and virtual reality. Princeton University Press., Princeton, N.J.
- Li and Mourikis, (2012) Li, M. and Mourikis, A. (2012). Improving the accuracy of EKF-based visual-inertial odometry. In Robotics and Automation (ICRA), 2012 IEEE International Conference on, pages 828–835.
- Li et al., (2014) Li, M., Yu, H., Zheng, X., and Mourikis, A. (2014). High-fidelity sensor modeling and self-calibration in vision-aided inertial navigation. In Robotics and Automation (ICRA), 2014 IEEE International Conference on, pages 409–416.
- Lupton and Sukkarieh, (2009) Lupton, T. and Sukkarieh, S. (2009). Efficient integration of inertial observations into visual slam without initialization. In IEEE/RSJ Int. Conf. on Intelligent Robots and Systems.
- Madyastha et al., (2011) Madyastha, V., Ravindra, V., Mallikarjunan, S., and Goyal, A. (2011). Extended Kalman Filter vs. Error State Kalman Filter for Aircraft Attitude Estimation. In AIAA Guidance, Navigation, and Control Conference, Portland, Oregon. American Institute of Aeronautics and Astronautics.
- Martinelli, (2012) Martinelli, A. (2012). Vision and imu data fusion: Closed-form solutions for attitude, speed, absolute scale, and bias determination. IEEE Transactions on Robotics, 28(1):44 –60.
- Piniés et al., (2007) Piniés, P., Lupton, T., Sukkarieh, S., and Tardós, J. D. (2007). Inertial aiding of inverse depth SLAM using a monocular camera. In Int. Conf. on Robotics and Automation.
- Roussillon et al., (2011) Roussillon, C., Gonzalez, A., Solà, J., Codol, J. M., Mansard, N., Lacroix, S., and Devy, M. (2011). RT-SLAM: A generic and real-time visual SLAM implementation. In Crowley, J., Draper, B., and Thonnat, M., editors, Computer Vision Systems, volume 6962 of Lecture Notes in Computer Science, pages 31–40. Springer Berlin Heidelberg.
- Shoemake, (1985) Shoemake, K. (1985). Animating rotation with quaternion curves. SIGGRAPH Comput. Graph., 19(3):245–254.
- Shuster, (1993) Shuster, M. D. (1993). Survey of attitude representations. Journal of the Astronautical Sciences, 41:439–517.
- Trawny and Roumeliotis, (2005) Trawny, N. and Roumeliotis, S. I. (2005). Indirect Kalman filter for 3D attitude estimation. Technical Report 2005-002, University of Minnesota, Dept. of Comp. Sci. & Eng.
- Yazell, (2009) Yazell, D. (2009). Origins of the unusual space shuttle quaternion definition. In Aerospace Sciences Meetings, pages –. American Institute of Aeronautics and Astronautics.