array(2) { ["docs"]=> array(10) { [0]=> array(10) { ["id"]=> string(3) "428" ["text"]=> string(77) "Visual Studio 2017 单独启动MSDN帮助(Microsoft Help Viewer)的方法" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(8) "DonetRen" ["tagsname"]=> string(55) "Visual Studio 2017|MSDN帮助|C#程序|.NET|Help Viewer" ["tagsid"]=> string(23) "[401,402,403,"300",404]" ["catesname"]=> string(0) "" ["catesid"]=> string(2) "[]" ["createtime"]=> string(10) "1511400964" ["_id"]=> string(3) "428" } [1]=> array(10) { ["id"]=> string(3) "427" ["text"]=> string(42) "npm -v;报错 cannot find module "wrapp"" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(4) "zzty" ["tagsname"]=> string(50) "node.js|npm|cannot find module "wrapp“|node" ["tagsid"]=> string(19) "[398,"239",399,400]" ["catesname"]=> string(0) "" ["catesid"]=> string(2) "[]" ["createtime"]=> string(10) "1511400760" ["_id"]=> string(3) "427" } [2]=> array(10) { ["id"]=> string(3) "426" ["text"]=> string(54) "说说css中pt、px、em、rem都扮演了什么角色" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(12) "zhengqiaoyin" ["tagsname"]=> string(0) "" ["tagsid"]=> string(2) "[]" ["catesname"]=> string(0) "" ["catesid"]=> string(2) "[]" ["createtime"]=> string(10) "1511400640" ["_id"]=> string(3) "426" } [3]=> array(10) { ["id"]=> string(3) "425" ["text"]=> string(83) "深入学习JS执行--创建执行上下文(变量对象,作用域链,this)" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(7) "Ry-yuan" ["tagsname"]=> string(33) "Javascript|Javascript执行过程" ["tagsid"]=> string(13) "["169","191"]" ["catesname"]=> string(0) "" ["catesid"]=> string(2) "[]" ["createtime"]=> string(10) "1511399901" ["_id"]=> string(3) "425" } [4]=> array(10) { ["id"]=> string(3) "424" ["text"]=> string(30) "C# 排序技术研究与对比" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(9) "vveiliang" ["tagsname"]=> string(0) "" ["tagsid"]=> string(2) "[]" ["catesname"]=> string(8) ".Net Dev" ["catesid"]=> string(5) "[199]" ["createtime"]=> string(10) "1511399150" ["_id"]=> string(3) "424" } [5]=> array(10) { ["id"]=> string(3) "423" ["text"]=> string(72) "【算法】小白的算法笔记:快速排序算法的编码和优化" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(9) "penghuwan" ["tagsname"]=> string(6) "算法" ["tagsid"]=> string(7) "["344"]" ["catesname"]=> string(0) "" ["catesid"]=> string(2) "[]" ["createtime"]=> string(10) "1511398109" ["_id"]=> string(3) "423" } [6]=> array(10) { ["id"]=> string(3) "422" ["text"]=> string(64) "JavaScript数据可视化编程学习(二)Flotr2,雷达图" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(7) "chengxs" ["tagsname"]=> string(28) "数据可视化|前端学习" ["tagsid"]=> string(9) "[396,397]" ["catesname"]=> string(18) "前端基本知识" ["catesid"]=> string(5) "[198]" ["createtime"]=> string(10) "1511397800" ["_id"]=> string(3) "422" } [7]=> array(10) { ["id"]=> string(3) "421" ["text"]=> string(36) "C#表达式目录树(Expression)" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(4) "wwym" ["tagsname"]=> string(0) "" ["tagsid"]=> string(2) "[]" ["catesname"]=> string(4) ".NET" ["catesid"]=> string(7) "["119"]" ["createtime"]=> string(10) "1511397474" ["_id"]=> string(3) "421" } [8]=> array(10) { ["id"]=> string(3) "420" ["text"]=> string(47) "数据结构 队列_队列实例:事件处理" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(7) "idreamo" ["tagsname"]=> string(40) "C语言|数据结构|队列|事件处理" ["tagsid"]=> string(23) "["246","247","248",395]" ["catesname"]=> string(12) "数据结构" ["catesid"]=> string(7) "["133"]" ["createtime"]=> string(10) "1511397279" ["_id"]=> string(3) "420" } [9]=> array(10) { ["id"]=> string(3) "419" ["text"]=> string(47) "久等了,博客园官方Android客户端发布" ["intro"]=> string(288) "目录 ECharts 异步加载 ECharts 数据可视化在过去几年中取得了巨大进展。开发人员对可视化产品的期望不再是简单的图表创建工具,而是在交互、性能、数据处理等方面有更高的要求。 chart.setOption({ color: [ " ["username"]=> string(3) "cmt" ["tagsname"]=> string(0) "" ["tagsid"]=> string(2) "[]" ["catesname"]=> string(0) "" ["catesid"]=> string(2) "[]" ["createtime"]=> string(10) "1511396549" ["_id"]=> string(3) "419" } } ["count"]=> int(200) } 222 视觉SLAM14讲笔记:ch2-4李群与李代数 - 爱码网

三维空间刚体运动

视觉SLAM14讲笔记:ch2-4李群与李代数

SLAM数学表达

xk=f(xk1,uk,wk)\boldsymbol{x}_{k}=f\left(\boldsymbol{x}_{k-1}, \boldsymbol{u}_{k}, \boldsymbol{w}_{k}\right)
这里uk\boldsymbol{u}_{k}是运动传感器的读数(有时也叫输入),wk\boldsymbol{w}_{k}为噪声。我们把它称为
运动方程

与运动方程相对应,还有一个观测方程。当小萝卜在xk\boldsymbol{x}_{k}位置上
看到某个路标点yj\boldsymbol{y}_{j},产生了一个观测数据zk,j\boldsymbol{z}_{k, j}
zk,j=h(yj,xk,vk,j)\boldsymbol{z}_{k, j}=h\left(\boldsymbol{y}_{j}, \boldsymbol{x}_{k}, \boldsymbol{v}_{k, j}\right)
这里vk,j\boldsymbol{v}_{k, j}是这次观测里的噪声。

旋转矩阵

点和向量,坐标系

三维空间中的某个向量的坐标可以用R3\mathbb{R}^{3}当中的三个数来描述。如果我们确定一个坐标系,那就可以谈论向量a 在这组基下的坐标了:
a=[e1,e2,e3][a1a2a3]=a1e1+a2e2+a3e3\boldsymbol{a}=\left[\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \boldsymbol{e}_{3}\right]\left[\begin{array}{l} a_{1} \\ a_{2} \\ a_{3} \end{array}\right]=a_{1} e_{1}+a_{2} e_{2}+a_{3} e_{3}

内积可以描述向量间的投影关系。而外积呢是这个样子:
a×b=[ijka1a2a3b1b2b3]=[a2b3a3b2a3b1a1b3a1b2a2b1]=[0a3a2a30a1a2a10]bab\boldsymbol{a} \times \boldsymbol{b}=\left[\begin{array}{ccc} \boldsymbol{i} & \boldsymbol{j} & \boldsymbol{k} \\ a_{1} & a_{2} & a_{3} \\ b_{1} & b_{2} & b_{3} \end{array}\right]=\left[\begin{array}{c} a_{2} b_{3}-a_{3} b_{2} \\ a_{3} b_{1}-a_{1} b_{3} \\ a_{1} b_{2}-a_{2} b_{1} \end{array}\right]=\left[\begin{array}{ccc} 0 & -a_{3} & a_{2} \\ a_{3} & 0 & -a_{1} \\ -a_{2} & a_{1} & 0 \end{array}\right] \boldsymbol{b} \triangleq \boldsymbol{a}^{\wedge} \boldsymbol{b}
我们引入了人符^号,把a\boldsymbol{a}写成一个矩阵。事实上是一个反对称矩阵 (Skew-symmetric),你可以将^记成一个反对称符号。这样就把外积a×ba \times b, 写成了矩阵与向量的乘法aba^{\wedge} b,把它变成了线性运算。

坐标系间的欧氏变换

相机运动是一个刚体运动,它保证了同一个向量在各个坐标系下的长度和夹角都不会 发生变化。这种变换称为欧氏变换。

首先来考虑旋转。我们设某个单位
正交基(e1,e2,e3)\left(e_{1}, e_{2}, e_{3}\right)经过一次旋转,变成了(e1,e2,e3)\left(e_{1}', e_{2}', e_{3}'\right)。那么,对于同一个向量a (注意该 向量并没有随着坐标系的旋转而发生运动),它在两个坐标系下的坐标为[a1,a2,a3]T\left[a_{1}, a_{2}, a_{3}\right]^{T}
[a1,a2,a3]T\left[a_{1}', a_{2}', a_{3}'\right]^{T}。根据坐标的定义,有:
[e1,e2,e3][a1a2a3]=[e1,e2,e3][a1a2a3]\left[e_{1}, e_{2}, e_{3}\right]\left[\begin{array}{c} a_{1} \\ a_{2} \\ a_{3} \end{array}\right]=\left[e_{1}^{\prime}, e_{2}^{\prime}, e_{3}^{\prime}\right]\left[\begin{array}{c} a_{1}^{\prime} \\ a_{2}^{\prime} \\ a_{3}^{\prime} \end{array}\right]
所以:
[a1a2a3]=[e1Te1e1Te2e1Te3e2Te1e2Te2e2Te3e3Te1e3Te2e3Te3][a1a2a3]Ra\left[\begin{array}{l} a_{1} \\ a_{2} \\ a_{3} \end{array}\right]=\left[\begin{array}{ccc} \boldsymbol{e}_{1}^{T} \boldsymbol{e}_{1}^{\prime} & \boldsymbol{e}_{1}^{T} \boldsymbol{e}_{2}^{\prime} & \boldsymbol{e}_{1}^{T} \boldsymbol{e}_{3}^{\prime} \\ \boldsymbol{e}_{2}^{T} \boldsymbol{e}_{1}^{\prime} & \boldsymbol{e}_{2}^{T} \boldsymbol{e}_{2}^{\prime} & \boldsymbol{e}_{2}^{T} \boldsymbol{e}_{3}^{\prime} \\ \boldsymbol{e}_{3}^{T} \boldsymbol{e}_{1}^{\prime} & \boldsymbol{e}_{3}^{T} \boldsymbol{e}_{2}^{\prime} & \boldsymbol{e}_{3}^{T} \boldsymbol{e}_{3}^{\prime} \end{array}\right]\left[\begin{array}{l} a_{1}^{\prime} \\ a_{2}^{\prime} \\ a_{3}^{\prime} \end{array}\right] \triangleq \boldsymbol{R} \boldsymbol{a}^{\prime}
可以说,矩阵R描述了旋转本身。因此它又称为旋转矩阵。

变换矩阵与齐次坐标

[a1]=[Rt0T1][a1]T[a1]\left[\begin{array}{l} \boldsymbol{a}^{\prime} \\ 1 \end{array}\right]=\left[\begin{array}{ll} \boldsymbol{R} & \boldsymbol{t} \\ \boldsymbol{0}^{T} & 1 \end{array}\right]\left[\begin{array}{l} \boldsymbol{a} \\ 1 \end{array}\right] \triangleq \boldsymbol{T}\left[\begin{array}{l} \boldsymbol{a} \\ 1 \end{array}\right]
我们把一个三维向量的末尾添加1,变成了四维向量,称为齐次
坐标。对于这个四维向量,我们可以把旋转和平移写在一个矩阵里面,使得整个关系变成
了线性关系。该式中,矩阵©称为变换矩阵(Transform Matrix)。我们暂时用N 表示a 的齐次坐标。

旋转向量和欧拉角

任意旋转都可以用一个旋转轴和一个旋转角来刻画。于是,我们可以使
用一个向量,其方向与旋转轴一致,而长度等于旋转角。这种向量,称为旋转向量(或轴角,Axis-Angle)。
而欧拉角则提供了一种非常直观的 方式来描述旋转—— 它使用了三个分离的转角,把一个旋转分解成三次绕不同轴的旋转。

  1. 绕物体的Z 轴旋转,得到偏航角yaw;
  2. 绕旋转之后的Y 轴旋转,得到俯仰角pitch;
  3. 绕旋转之后的X 轴旋转,得到滚转角roll。

四元数

旋转矩阵用九个量描述三自由度的旋转,具有冗余性;欧拉角和旋转向量是紧凑的, 但具有奇异性。事实上,我们找不到不带奇异性的三维向量描述方式。一个四元数q 拥有一个实部和三个虚部。
q=q0+q1i+q2j+q3k\boldsymbol{q}=q_{0}+q_{1} i+q_{2} j+q_{3} k
其中i,j,k为四元数的三个虚部。这三个虚部满足关系式:
{i2=j2=k2=1ij=k,ji=kjk=i,kj=iki=j,ik=j\left\{\begin{array}{l} i^{2}=j^{2}=k^{2}=-1 \\ i j=k, j i=-k \\ j k=i, k j=-i \\ k i=j, i k=-j \end{array}\right.
我们知道单位四元数能够表达三维空间的旋转。这种表达方式和旋转矩
阵、旋转向量有什么关系呢?
q=[cosθ2,nxsinθ2,nysinθ2,nzsinθ2]T\boldsymbol{q}=\left[\cos \frac{\theta}{2}, n_{x} \sin \frac{\theta}{2}, n_{y} \sin \frac{\theta}{2}, n_{z} \sin \frac{\theta}{2}\right]^{T}
反之,我们亦可从单位四元数中计算出对应旋转轴与夹角:
{θ=2arccosq0[nx,ny,nz]T=[q1,q2,q3]T/sinθ2\left\{\begin{array}{l} \theta=2 \arccos q_{0} \\ {\left[n_{x}, n_{y}, n_{z}\right]^{T}=\left[q_{1}, q_{2}, q_{3}\right]^{T} / \sin \frac{\theta}{2}} \end{array}\right.

李群与李代数

李群与李代数的定义

上一讲,我们介绍了三维世界中刚体运动的描述方式,包括旋转矩阵、旋转向量、欧 拉角、四元数等若干种方式。我们重点介绍了旋转的表示,但是在SLAM中,除了表示之外,我们还要对它们进行估计和优化。因为在SLAM中位姿是未知的,而我们需要解决什么样的相机位姿最符合当前观测数据这样的问题。一种典型的方式是把它构建成一个优化问题。

(Group)是一种集合加上一种运算的代数结构。我们把集合记作4 , 运算记作・, 那么群可以记作G=(A,)G=(A, \cdot), 群要求这个运算满足以下几个条件:

  1. 封闭性: a1,a2A,a1a2A\quad \forall a_{1}, a_{2} \in A, \quad a_{1} \cdot a_{2} \in A

  2. 结合律: a1,a2,a3A,(a1a2)a3=a1(a2a3)\quad \forall a_{1}, a_{2}, a_{3} \in A, \quad\left(a_{1} \cdot a_{2}\right) \cdot a_{3}=a_{1} \cdot\left(a_{2} \cdot a_{3}\right)

  3. 幺元: a0A,s.t.aA,a0a=aa0=a\quad \exists a_{0} \in A, \quad s.t.\quad \forall a \in A, \quad a_{0} \cdot a=a \cdot a_{0}=a

  4. 逆: aA,a1A,\quad \forall a \in A, \quad \exists a^{-1} \in A, \quad s.t. aa1=a0\quad a \cdot a^{-1}=a_{0}

每个李群都有与之对应的李代数。李代数描述了李群的局部性质。通用的李代数的定义如下:
李代数由一个集合V\mathbb{V} , —个数域F\mathbb{F} 和一个二元运算[,]组成。如果它们满足以下几条
性质,称(V,F,[,])(\mathbb{V}, \mathbb{F},[,])为一个李代数,记作g\mathfrak{g}

  1. 封闭性 X,YV,[X,Y]V\quad \forall \boldsymbol{X}, \boldsymbol{Y} \in \mathbb{V},[\boldsymbol{X}, \boldsymbol{Y}] \in \mathbb{V}
  2. 双线性 X,Y,ZV,a,bF,\quad \forall \boldsymbol{X}, \boldsymbol{Y}, \boldsymbol{Z} \in \mathbb{V}, a, b \in \mathbb{F}, 有 :
    $
    [a \boldsymbol{X}+b \boldsymbol{Y}, \boldsymbol{Z}]=a[\boldsymbol{X}, \boldsymbol{Z}]+b[\boldsymbol{Y}, \boldsymbol{Z}], \quad[\boldsymbol{Z}, a \boldsymbol{X}+b \boldsymbol{Y}]=a[\boldsymbol{Z}, \boldsymbol{X}]+b[\boldsymbol{Z}, \boldsymbol{Y}]
    $
  3. 自反性 XV,[X,X]=0\quad \forall \boldsymbol{X} \in \mathbb{V},[\boldsymbol{X}, \boldsymbol{X}]=\mathbf{0}
  4. 雅可比等价 X,Y,ZV,[X,[Y,Z]]+[Z,[Y,X]]+[Y,[Z,X]]=0\quad \forall \boldsymbol{X}, \boldsymbol{Y}, \boldsymbol{Z} \in \mathbb{V},[\boldsymbol{X},[\boldsymbol{Y}, \boldsymbol{Z}]]+[\boldsymbol{Z},[\boldsymbol{Y}, \boldsymbol{X}]]+[\boldsymbol{Y},[\boldsymbol{Z}, \boldsymbol{X}]]=\mathbf{0}

其中二元运算被称为李括号。从表面上来看,李代数所需要的性质还是挺多的。相比于群中的较为简单的二元运算,李括号表达了两个元素的差异。它不要求结合律,而要求元素和自己做李括号之后为零的性质。
作为例子,三维向量上定义的叉积是一种李括号,g=(R3,R,×)\mathfrak{g}=\left(\mathbb{R}^{3}, \mathbb{R}, \times\right)构成了一个李代数。

指数与对数映射

SO(3) 上 的 指 数 映 射

根据前面的推导,每个 ϕ\phi 都可以生成一个反对称矩阵:
Φ=ϕ=[0ϕ3ϕ2ϕ30ϕ1ϕ2ϕ10]R3×3\mathbf{\Phi}=\boldsymbol{\phi}^{\wedge}=\left[\begin{array}{ccc} 0 & -\phi_{3} & \phi_{2} \\ \phi_{3} & 0 & -\phi_{1} \\ -\phi_{2} & \phi_{1} & 0 \end{array}\right] \in \mathbb{R}^{3 \times 3}
在此定义下,两个向量 ϕ1,ϕ2\phi_{1}, \phi_{2} 的李括号为:
[ϕ1,ϕ2]=(Φ1Φ2Φ2Φ1) \left[\phi_{1}, \phi_{2}\right]=\left(\Phi_{1} \Phi_{2}-\Phi_{2} \Phi_{1}\right)^{\vee}
李群为:
so(3)={ϕR3,Φ=ϕR3×3}\mathfrak{s o}(3)=\left\{\boldsymbol{\phi} \in \mathbb{R}^{3}, \boldsymbol{\Phi}=\boldsymbol{\phi}^{\wedge} \in \mathbb{R}^{3 \times 3}\right\}

现在,我们说,R 是某个相机的旋转,它会随时间连续地变化,即为时间的函数:R(t)由于它仍是旋转矩阵,有
R(t)R(t)T=I\boldsymbol{R}(t) \boldsymbol{R}(t)^{T}=\boldsymbol{I}
在等式两边对时间求导:
R˙(t)R(t)T=(R˙(t)R(t)T)T\dot{\boldsymbol{R}}(t) \boldsymbol{R}(t)^{T}=-\left(\dot{\boldsymbol{R}}(t) \boldsymbol{R}(t)^{T}\right)^{T}
引入了^符号,将一个向量变成了反对称矩阵。同理,对于任意反对称矩阵,我们亦能找 到一个与之对应的向量。把这个运算用符号^\vee表示:
a=A=[0a3a2a30a1a2a10],A=a\boldsymbol{a}^{\wedge}=\boldsymbol{A}=\left[\begin{array}{ccc} 0 & -a_{3} & a_{2} \\ a_{3} & 0 & -a_{1} \\ -a_{2} & a_{1} & 0 \end{array}\right], \quad \boldsymbol{A}^{\vee}=\boldsymbol{a}
于是,由于 R˙(t)R(t)T\dot{\boldsymbol{R}}(t) \boldsymbol{R}(t)^{T} 是一个反对称矩阵,我们可以找到一个三维向量 ϕ(t)R3\boldsymbol{\phi}(t) \in \mathbb{R}^{3}
之对应。于是有:
R˙(t)R(t)T=ϕ(t)\dot{\boldsymbol{R}}(t) \boldsymbol{R}(t)^{T}=\boldsymbol{\phi}(t)^{\wedge}
等式两边右乘 R(t),\boldsymbol{R}(t), 由于 R\boldsymbol{R} 为正交阵,有
R˙(t)=ϕ(t)R(t)=[0ϕ3ϕ2ϕ30ϕ1ϕ2ϕ10]R(t)\dot{\boldsymbol{R}}(t)=\phi(t)^{\wedge} \boldsymbol{R}(t)=\left[\begin{array}{ccc} 0 & -\phi_{3} & \phi_{2} \\ \phi_{3} & 0 & -\phi_{1} \\ -\phi_{2} & \phi_{1} & 0 \end{array}\right] \boldsymbol{R}(t)
上式是一个关于 R 的微分方程,而且我们知道初始值 R(0)=I,\boldsymbol{R}(0)=\boldsymbol{I}, 解之,得:
R(t)=exp(ϕ0t)\boldsymbol{R}(t)=\exp \left(\phi_{0}^{\wedge} t\right)
exp(ϕ)\exp \left(\phi^{\wedge}\right) 是如何计算的?它是一个矩阵的指数,在李群和李 代数中,称为指数映射(ExponentialMap)。
任意矩阵的指数映射可以写成一个泰勒展开,但是只有在收敛的情况下才会有结果, 其结果仍是一个矩阵。
exp(ϕ)=n=01n!(ϕ)n\exp \left(\phi^{\wedge}\right)=\sum_{n=0}^{\infty} \frac{1}{n !}\left(\phi^{\wedge}\right)^{n}
我们来仔细推导一下这个定义。由于 ϕ\phi 是三维向量,我们可以定义它的模长和它的
方向,分别记作 θ\thetaa,\boldsymbol{a}, 于是有 ϕ=θa\boldsymbol{\phi}=\theta \boldsymbol{a} \circ 这里 a\boldsymbol{a} 是一个长度为 1 的方向向量。首先,对
a,\boldsymbol{a}^{\wedge}, 有以下两条性质:
aa=aaTIaaa=a\begin{aligned} &\boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}=\boldsymbol{a} \boldsymbol{a}^{T}-\boldsymbol{I}\\ &\boldsymbol{a ^ { \wedge }} \boldsymbol{a ^ { \wedge }} \boldsymbol{a}^{\wedge}=-\boldsymbol{a}^{\wedge} \end{aligned}
它们提供了处理a\boldsymbol{a}^{\wedge}高阶项的方法。利用这两个性质, 我们可以把指数映射写成:
exp(ϕ)=exp(θa)=n=01n!(θa)n=I+θa+12!θ2aa+13!θ3aaa+14!θ4(a)4+=aaTaa+θa+12!θ2aa13!θ3a14!θ4(a)2+=aaT+(θ13!θ3+15!θ5)a(112!θ2+14!θ4)aa=aa+I+sinθacosθaa=(1cosθ)aa+I+sinθa=cosθI+(1cosθ)aaT+sinθa\begin{aligned} \exp \left(\phi^{\wedge}\right) &=\exp \left(\theta \boldsymbol{a}^{\wedge}\right)=\sum_{n=0}^{\infty} \frac{1}{n !}\left(\theta \boldsymbol{a}^{\wedge}\right)^{n} \\ &=\boldsymbol{I}+\theta \boldsymbol{a}^{\wedge}+\frac{1}{2 !} \theta^{2} \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}+\frac{1}{3 !} \theta^{3} \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}+\frac{1}{4 !} \theta^{4}\left(\boldsymbol{a}^{\wedge}\right)^{4}+\ldots \\ &=\boldsymbol{a} \boldsymbol{a}^{T}-\boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}+\theta \boldsymbol{a}^{\wedge}+\frac{1}{2 !} \theta^{2} \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}-\frac{1}{3 !} \theta^{3} \boldsymbol{a}^{\wedge}-\frac{1}{4 !} \theta^{4}\left(\boldsymbol{a}^{\wedge}\right)^{2}+\ldots \\ &=\boldsymbol{a} \boldsymbol{a}^{T}+\left(\theta-\frac{1}{3 !} \theta^{3}+\frac{1}{5 !} \theta^{5}-\ldots\right) \boldsymbol{a}^{\wedge}-\left(1-\frac{1}{2 !} \theta^{2}+\frac{1}{4 !} \theta^{4}-\ldots\right) \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge} \\ &=\boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}+\boldsymbol{I}+\sin \theta \boldsymbol{a}^{\wedge}-\cos \theta \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge} \\ &=(1-\cos \theta) \boldsymbol{a}^{\wedge} \boldsymbol{a}^{\wedge}+\boldsymbol{I}+\sin \theta \boldsymbol{a}^{\wedge} \\ &=\cos \theta \boldsymbol{I}+(1-\cos \theta) \boldsymbol{a} \boldsymbol{a}^{T}+\sin \theta \boldsymbol{a}^{\wedge} \end{aligned}
最后我们得到罗德里格斯公式:
exp(θa)=cosθI+(1cosθ)aaT+sinθa\exp \left(\theta \boldsymbol{a}^{\wedge}\right)=\cos \theta \boldsymbol{I}+(1-\cos \theta) \boldsymbol{a} \boldsymbol{a}^{T}+\sin \theta \boldsymbol{a}^{\wedge}

SE (3)上的指数映射

对于 SE(3),S E(3), 它也有对应的李代数 se(3)\mathfrak{se}(3) \circ 为省略篇幅,我们就不描述如何引出 se(3)\mathfrak{se}(3)
了。与 so(3)\mathfrak{so}(3) 相似,se(3)\mathfrak{se}(3) 位于 R6\mathbb{R}^{6} 空间中:
se(3)={ξ=[ρϕ]R6,ρR3,ϕso(3),ξ=[ϕρ0T0]R4×4}\mathfrak{se}(3)=\left\{\boldsymbol{\xi}=\left[\begin{array}{c}\boldsymbol{\rho} \\ \boldsymbol{\phi}\end{array}\right] \in \mathbb{R}^{6}, \boldsymbol{\rho} \in \mathbb{R}^{3}, \boldsymbol{\phi} \in \mathfrak{s o}(3), \boldsymbol{\xi}^{\wedge}=\left[\begin{array}{cc}\boldsymbol{\phi}^{\wedge} & \boldsymbol{\rho} \\ \mathbf{0}^{T} & 0\end{array}\right] \in \mathbb{R}^{4 \times 4}\right\}
我们把每个 se(3)\mathfrak{se}(3) 元素记作 ξ,\boldsymbol{\xi}, 它是一个六维向量。前三维为平移,记作 ρ;\boldsymbol{\rho} ; 后三维为旋转,记作 ϕ,\phi, 实质上是 so(3)\mathfrak{so}(3) 元素。同时,我们拓展了^\wedge符号的含义。在 se(3)\mathfrak{se}(3) 中, 同样使用 ^ 符号,将一个六维向量转换成四维矩阵,但这里不再表示反对称:
ξ=[ϕρ0T0]R4×4 \boldsymbol{\xi}^{\wedge}=\left[\begin{array}{cc} \boldsymbol{\phi}^{\wedge} & \boldsymbol{\rho} \\ \mathbf{0}^{T} & 0 \end{array}\right] \in \mathbb{R}^{4 \times 4}
下面我们来介绍 se(3)\mathfrak{se}(3) 上的指数映射。为了节省篇幅,我们不再像 so(3)\mathfrak{so}(3) 那样详细推 导指数映射。se(3)\mathfrak{se}(3) 上的指数映射形式如下:
exp(ξ)=[n=01n!(ϕ)nn=01(n+1)!(ϕ)nρ0T1][RJρ0T1]=T\begin{aligned} \exp \left(\boldsymbol{\xi}^{\wedge}\right) &=\left[\begin{array}{ccc} \sum_{n=0}^{\infty} \frac{1}{n !}\left(\phi^{\wedge}\right)^{n} & \sum_{n=0}^{\infty} \frac{1}{(n+1) !}\left(\phi^{\wedge}\right)^{n} \boldsymbol{\rho} \\ \mathbf{0}^{T} & 1 \end{array}\right] \\ & \triangleq\left[\begin{array}{cc} \boldsymbol{R} & \boldsymbol{J} \boldsymbol{\rho} \\ \mathbf{0}^{T} & 1 \end{array}\right]=\boldsymbol{T} \end{aligned}
如果你有耐心,可以照着 so(3)\mathfrak{so}(3) 上的做法推导,把 exp 进行泰勒展开推导此式。从结 果上看,ξ\xi的指数映射左上角的 R 是我们熟知的 so(3)\mathfrak{so}(3) 中的元素,与 se(3)\mathfrak{se}(3) 当中的旋转部
ϕ\phi 对应。而右上角的 J 则可整理为(设 ϕ=θa\phi=\theta a):
J=sinθθI+(1sinθθ)aaT+1cosθθaJ=\frac{\sin \theta}{\theta} I+\left(1-\frac{\sin \theta}{\theta}\right) a a^{T}+\frac{1-\cos \theta}{\theta} a^{\wedge}
视觉SLAM14讲笔记:ch2-4李群与李代数

李代数求导与扰动模型

BCH公式与近似形式

BCH公式告诉我们,当处理两个矩阵指数之积时,它们会产生一些由
李括号组成的余项。
In the particular case of SO(3)S O(3), we can show that
ln(C1C2)=ln(exp(ϕ1)exp(ϕ2))=ϕ1+ϕ2+12ϕ1ϕ2+112ϕ1ϕ1ϕ2+112ϕ2ϕ2ϕ1+ \begin{array}{l} \ln \left(\mathbf{C}_{1} \mathbf{C}_{2}\right)^{\vee}=\ln \left(\exp \left(\phi_{1}^{\wedge}\right) \exp \left(\phi_{2}^{\wedge}\right)\right)^{\vee} \\ \quad=\phi_{1}+\phi_{2}+\frac{1}{2} \phi_{1}^{\wedge} \phi_{2}+\frac{1}{12} \phi_{1}^{\wedge} \phi_{1}^{\wedge} \phi_{2}+\frac{1}{12} \phi_{2}^{\wedge} \phi_{2}^{\wedge} \phi_{1}+\cdots \end{array}
where C1=exp(ϕ1),C2=exp(ϕ2)SO(3).\mathbf{C}_{1}=\exp \left(\phi_{1}^{\wedge}\right), \mathbf{C}_{2}=\exp \left(\phi_{2}^{\wedge}\right) \in S O(3) . Alternatively, if we as-
sume that ϕ1\phi_{1} or ϕ2\phi_{2} is small, then using the approximate BCH formulas we can show that
ln(C1C2)=ln(exp(ϕ1)exp(ϕ2)){J(ϕ2)1ϕ1+ϕ2 if ϕ1 small ϕ1+Jr(ϕ1)1ϕ2 if ϕ2 small  \begin{aligned} \ln \left(\mathbf{C}_{1} \mathbf{C}_{2}\right)^{\vee}=\ln \left(\exp \left(\phi_{1}^{\wedge}\right) \exp \left(\phi_{2}^{\wedge}\right)\right)^{\vee} \approx\left\{\begin{aligned} \mathbf{J}_{\ell}\left(\phi_{2}\right)^{-1} \phi_{1}+\phi_{2} & \text { if } \phi_{1} \text { small } \\ \phi_{1}+\mathbf{J}_{r}\left(\phi_{1}\right)^{-1} \phi_{2} & \text { if } \phi_{2} \text { small } \end{aligned}\right. \end{aligned}
本书以左乘为例。左乘 BCH 近似雅可比 JlJ_{l}实际上就是SE(3)SE (3)上的指数映射中的JJ
Jl=J=sinθθI+(1sinθθ)aaT+1cosθθa \boldsymbol{J}_{l}=\boldsymbol{J}=\frac{\sin \theta}{\theta} \boldsymbol{I}+\left(1-\frac{\sin \theta}{\theta}\right) \boldsymbol{a} \boldsymbol{a}^{T}+\frac{1-\cos \theta}{\theta} \boldsymbol{a}^{\wedge}
它的逆为:
Jl1=θ2cotθ2I+(1θ2cotθ2)aaTθ2a \boldsymbol{J}_{l}^{-1}=\frac{\theta}{2} \cot \frac{\theta}{2} \boldsymbol{I}+\left(1-\frac{\theta}{2} \cot \frac{\theta}{2}\right) \boldsymbol{a} \boldsymbol{a}^{T}-\frac{\theta}{2} \boldsymbol{a}^{\wedge}
而右乘雅可比仅需要对自变量取负号即可:
Jr(ϕ)=Jl(ϕ) \boldsymbol{J}_{r}(\boldsymbol{\phi})=\boldsymbol{J}_{l}(-\boldsymbol{\phi})

假定对某个旋转 R,对应的李代数为 ϕ\phi 。我们给它左乘一个微小旋转,记作 ΔR\Delta \boldsymbol{R} 对应的李代数为 Δϕ\Delta \phi 。 那么,在李群上,得到的结果就是 ΔRR,\Delta \boldsymbol{R} \cdot \boldsymbol{R}, 而在李代数上,根据 BCH 近似,为: Jl1(ϕ)Δϕ+ϕJ_{l}^{-1}(\phi) \Delta \phi+\phi 。合并起来,可以简单地写成:
exp(Δϕ)exp(ϕ)=exp((ϕ+Jl1(ϕ)Δϕ)) \exp \left(\Delta \phi^{\wedge}\right) \exp \left(\phi^{\wedge}\right)=\exp \left(\left(\phi+J_{l}^{-1}(\phi) \Delta \phi\right)^{\wedge}\right)
反之,如果我们在李代数上进行加法,让一个 ϕ\phi 加上 Δϕ,\Delta \phi, 那么可以近似为李群上带 左右雅可比的乘法:
exp((ϕ+Δϕ))=exp((JlΔϕ))exp(ϕ)=exp(ϕ)exp((JrΔϕ)) \exp \left((\phi+\Delta \phi)^{\wedge}\right)=\exp \left(\left(J_{l} \Delta \phi\right)^{\wedge}\right) \exp \left(\phi^{\wedge}\right)=\exp \left(\phi^{\wedge}\right) \exp \left(\left(J_{r} \Delta \phi\right)^{\wedge}\right)
这将为之后李代数上的做微积分提供了理论基础。同样的,对于 SE(3),亦有类似的 BCH 近似公式:
exp(Δξ)exp(ξ)exp((Jl1Δξ+ξ)) \exp \left(\Delta \xi^{\wedge}\right) \exp \left(\boldsymbol{\xi}^{\wedge}\right) \approx \exp \left(\left(\mathcal{J}_{l}^{-1} \Delta \boldsymbol{\xi}+\boldsymbol{\xi}\right)^{\wedge}\right)
exp(ξ)exp(Δξ)exp((Jr1Δξ+ξ)) \exp \left(\boldsymbol{\xi}^{\wedge}\right) \exp \left(\Delta \boldsymbol{\xi}^{\wedge}\right) \approx \exp \left(\left(\mathcal{J}_{r}^{-1} \Delta \boldsymbol{\xi}+\boldsymbol{\xi}\right)^{\wedge}\right)

SO(3)李代数上的求导

下面我们来讨论一个带有李代数的函数,如何关于该李代数求导的问题。该问题有很
强的实际背景。在 SLAM中,我们要估计一个相机的位置和姿态,该位姿是由SO(3)SO(3)的旋转矩阵或SE(3)SE(3)上的变换矩阵描述的。不妨设某个时刻小萝卜的位姿为它观察到了一个世界坐标位于pp的点,产生了一个观测数据zz 。那么,由坐标变换关系知:
z=Tp+wz=T p+w
其中ww为噪声。我们通常会计算理想的观测与实际数据的误差:
e=zTp e=z-T p
假设一共有 N 个这样的路标点和观测,于是就有 NN 个上式。那么,对小萝卜的位姿估计,相当于是寻找一个最优的 TT,使得整体误差最小化:
minTJ(T)=i=1NziTpi22 \min _{\boldsymbol{T}} J(\boldsymbol{T})=\sum_{i=1}^{N}\left\|\boldsymbol{z}_{i}-\boldsymbol{T} \boldsymbol{p}_{i}\right\|_{2}^{2}
我们经常会构建与位姿有关的函数,然后讨论该函数关于位姿的导数,以调整当前的估计值。然而,SO(3)SO(3),SE(3)SE(3)上并没有良好定义的加法,它们只是群。如果我们把TT当成一个普通矩阵来处理优化,那就必须对它加以约束。而从李代数角度来说,由于李代数由向量组成,具有良好的加法运算。因此,使用李代数解决求导问题的思路分为两种:

  1. 用李代数表示姿态,然后对根据李代数加法来对李代数求导。
  2. 对李群左乘或右乘微小扰动,然后对该扰动求导,称为左扰动和右扰动模型。

首先,考虑SO(3)SO(3)上的情况。假设我们对一个空间点pp进行了旋转,得到了RpRp。现在,要计算旋转之后点的坐标相对于旋转的导数,由于 SO(3) 没有加法,所以该导数无法按?导数的定义进行计算。设 R 对应的李代数为 ϕ,\boldsymbol{\phi}, 我们转而计算:
(exp(ϕ)p)ϕ \frac{\partial\left(\exp \left(\phi^{\wedge}\right) p\right)}{\partial \phi}
按照导数的定义,有:
(exp(ϕ)p)ϕ=limδϕ0exp((ϕ+δϕ))pexp(ϕ)pδϕ=limδϕ0exp((Jlδϕ))exp(ϕ)pexp(ϕ)pδϕ=limδϕ0(I+(Jlδϕ))exp(ϕ)pexp(ϕ)pδϕ=limδϕ0(Jlδϕ)exp(ϕ)pδϕ=limδϕ0(exp(ϕ)p)Jlδϕδϕ=(Rp)Jl \begin{aligned} \frac{\partial\left(\exp \left(\phi^{\wedge}\right) \boldsymbol{p}\right)}{\partial \phi} &=\lim _{\delta \phi \rightarrow 0} \frac{\exp \left((\phi+\delta \phi)^{\wedge}\right) \boldsymbol{p}-\exp \left(\phi^{\wedge}\right) \boldsymbol{p}}{\delta \phi} \\ &=\lim _{\delta \phi \rightarrow 0} \frac{\exp \left(\left(\boldsymbol{J}_{l} \delta \phi\right)^{\wedge}\right) \exp \left(\boldsymbol{\phi}^{\wedge}\right) \boldsymbol{p}-\exp \left(\phi^{\wedge}\right) \boldsymbol{p}}{\delta \boldsymbol{\phi}} \\ &=\lim _{\delta \phi \rightarrow 0} \frac{\left(\boldsymbol{I}+\left(\boldsymbol{J}_{l} \delta \boldsymbol{\phi}\right)^{\wedge}\right) \exp \left(\boldsymbol{\phi}^{\wedge}\right) \boldsymbol{p}-\exp \left(\phi^{\wedge}\right) \boldsymbol{p}}{\delta \phi} \\ &=\lim _{\delta \phi \rightarrow 0} \frac{\left(\boldsymbol{J}_{l} \delta \boldsymbol{\phi}\right)^{\wedge} \exp \left(\phi^{\wedge}\right) \boldsymbol{p}}{\delta \boldsymbol{\phi}} \\ &=\lim _{\delta \phi \rightarrow 0} \frac{-\left(\exp \left(\phi^{\wedge}\right) \boldsymbol{p}\right)^{\wedge} \boldsymbol{J}_{l} \delta \boldsymbol{\phi}}{\delta \boldsymbol{\phi}}=-(\boldsymbol{R} \boldsymbol{p})^{\wedge} \boldsymbol{J}_{l} \end{aligned}
第二行的近似为BCH线性近似,第三行为泰勒展开舍去高阶项后近似,第四行至第五行将反对称符号看作叉积,交换之后变号。于是,我们推导了旋转后的点相对于李代数的导数:
(Rp)ϕ=(Rp)Jl\frac{\partial(\boldsymbol{R} \boldsymbol{p})}{\partial \boldsymbol{\phi}}=(-\boldsymbol{R} \boldsymbol{p})^{\wedge} \boldsymbol{J}_{l}

扰动模型(左乘)

另一种求导方式,是对$ R $进行一次扰动 ΔR\Delta \boldsymbol{R}。 这个扰动可以乘在左边也可以乘在右边,最后结果会有一点儿微小的差异,我们以左扰动为例。设左扰动 ΔR\Delta \boldsymbol{R} 对应的李代数 为φ\varphi。然后,对 φ\varphi 求导,即:
(Rp)φ=limφ0exp(φ)exp(ϕ)pexp(ϕ)pφlimφ0(1+φ)exp(ϕ)pexp(ϕ)pφ=limφ0φRpφ=limφ0(Rp)φφ=(Rp) \begin{aligned} \frac{\partial(\boldsymbol{R} \boldsymbol{p})}{\partial \boldsymbol{\varphi}} &=\lim _{\boldsymbol{\varphi} \rightarrow 0} \frac{\exp \left(\boldsymbol{\varphi}^{\wedge}\right) \exp \left(\boldsymbol{\phi}^{\wedge}\right) \boldsymbol{p}-\exp \left(\boldsymbol{\phi}^{\wedge}\right) \boldsymbol{p}}{\varphi} \\ & \approx \lim _{\boldsymbol{\varphi} \rightarrow 0} \frac{\left(1+\boldsymbol{\varphi}^{\wedge}\right) \exp \left(\boldsymbol{\phi}^{\wedge}\right) \boldsymbol{p}-\exp \left(\boldsymbol{\phi}^{\wedge}\right) \boldsymbol{p}}{\varphi} \\ &=\lim _{\boldsymbol{\varphi} \rightarrow 0} \frac{\boldsymbol{\varphi}^{\wedge} \boldsymbol{R} \boldsymbol{p}}{\boldsymbol{\varphi}}=\lim _{\varphi \rightarrow 0} \frac{-(\boldsymbol{R} \boldsymbol{p})^{\wedge} \boldsymbol{\varphi}}{\boldsymbol{\varphi}}=-(\boldsymbol{R} \boldsymbol{p})^{\wedge} \end{aligned}
可见,扰动模型相比于直接对李代数求导,省去了一个雅可比JlJ_l的计算。这使得扰动模型更为实用。请读者务必理解这里的求导运算,这在位姿估计当中具有重要的意义。

SE(3)上的李代数求导

最后,我们给出 SE(3) 上的扰动模型,而直接李代数上的求导就不再介绍了。假设 恭空间点 p\boldsymbol{p} 经过一次变换 T\boldsymbol{T} (对应李代数为 ξ\boldsymbol{\xi}), 得到 Tp\boldsymbol{T} \boldsymbol{p} 。现在,给 T\boldsymbol{T} 左乘一个扰动 ΔT=exp(δξ),\Delta \boldsymbol{T}=\exp \left(\delta \boldsymbol{\xi}^{\wedge}\right), 我们设扰动项的李代数为 δξ=[δρ,δϕ]T,\delta \boldsymbol{\xi}=[\delta \boldsymbol{\rho}, \delta \boldsymbol{\phi}]^{T}, \quad 那么:
(Tp)δξ=limδξ0exp(δξ)exp(ξ)pexp(ξ)pδξlimδξ0(I+δξ)exp(ξ)pexp(ξ)pδξ=limδξ0δξexp(ξ)pδξ=limδξ0[δϕδρ0T0][Rp+t1]δξ\begin{aligned} \frac{\partial(\boldsymbol{T} \boldsymbol{p})}{\partial \delta \boldsymbol{\xi}} &=\lim _{\delta \boldsymbol{\xi} \rightarrow 0} \frac{\exp \left(\delta \boldsymbol{\xi}^{\wedge}\right) \exp \left(\boldsymbol{\xi}^{\wedge}\right) \boldsymbol{p}-\exp \left(\boldsymbol{\xi}^{\wedge}\right) \boldsymbol{p}}{\delta \boldsymbol{\xi}} \\ & \approx \lim _{\delta \boldsymbol{\xi} \rightarrow 0} \frac{\left(\boldsymbol{I}+\delta \boldsymbol{\xi}^{\wedge}\right) \exp \left(\boldsymbol{\xi}^{\wedge}\right) \boldsymbol{p}-\exp \left(\boldsymbol{\xi}^{\wedge}\right) \boldsymbol{p}}{\delta \boldsymbol{\xi}} \\ &=\lim _{\delta \boldsymbol{\xi} \rightarrow 0} \frac{\delta \boldsymbol{\xi}^{\wedge} \exp \left(\boldsymbol{\xi}^{\wedge}\right) \boldsymbol{p}}{\delta \boldsymbol{\xi}}\\ &=\lim _{\delta \boldsymbol{\xi} \rightarrow 0} \frac{\left[\begin{array}{cc} \delta \boldsymbol{\phi}^{\wedge} & \delta \boldsymbol{\rho} \\ \mathbf{0}^{T} & 0 \end{array}\right]\left[\begin{array}{c} \boldsymbol{R} \boldsymbol{p}+\boldsymbol{t} \\ 1 \end{array}\right]}{\delta \boldsymbol{\xi}} \end{aligned}

=limδξ0[δϕ(Rp+t)+δρ0]δξ=[I(Rp+t)0T0T](Tp)=\lim _{\delta \boldsymbol{\xi} \rightarrow 0} \frac{\left[\begin{array}{c} \delta \phi^{\wedge}(\boldsymbol{R} \boldsymbol{p}+\boldsymbol{t})+\delta \boldsymbol{\rho} \\ 0 \end{array}\right]}{\delta \boldsymbol{\xi}}=\left[\begin{array}{cc} \boldsymbol{I} & -(\boldsymbol{R} \boldsymbol{p}+\boldsymbol{t})^{\wedge} \\ \mathbf{0}^{T} & \mathbf{0}^{T} \end{array}\right] \triangleq(\boldsymbol{T} \boldsymbol{p})^{\odot}

相关文章:

  • 2019-03-13
  • 2020-04-30
  • 2019-03-23
  • 2021-08-09
  • 2019-03-29
  • 2019-07-24
  • 2021-09-27
  • 2018-12-17
猜你喜欢
  • 2021-12-20
  • 2021-11-14
  • 2019-11-30
  • 2021-12-12
  • 2021-08-04
  • 2019-11-05
  • 2020-12-27
  • 2021-12-28
相关资源
相似解决方案