MathSword教程7.启发式优化求解器

启发式优化算法有比较强的全局寻优能力,因此对于非线性优化(非线性拟合、非线性方程组求解等等)问题可以优先尝试使用启发式优化求解器。

目前MathSword内置自研SC启发式优化算法、以及基于龙格库塔法的RUN[1]算法,这两种算法是我目前测试过智能算法中效果最好的两种算法。这里将启发式优化求解的使用方法演示如下:

 

参考

[1] Chen H , Ahmadianfar I , Heidari A A , et al. RUN beyond the metaphor: An efficient optimization algorithm based on Runge Kutta method[J]. Expert Systems with Application, 2021(181-Nov.).

MathSword教程8.求解线性以及非线性规划问题

1 说明

MathSword提供一些特定函数或者窗口来求解线性、非线性方程组,或者线性、非线性规划问题。对于自由度比较高的非线性问题求解,会涉及到程序语言编写问题。为了使得求解的问题,在软件输入上更加接近数学语言,因此MathSword专门提供一个规划求解窗口来按指定规则书写对应问题,即可以求解。对于旅行商问题、指派问题、背包问题、最短路径问题,由于不方便写成对应表达式,建议使用专门的求解函数进行求解。

2 教程

2.1 求解线性方程

求解如下线性方程组,其中\(x_1,x_2,x_3\)为浮点变量

$$\begin{cases}
-49x_1+82x_2+36x_3&=14.67
\\
92x_1+6x_2+x_3&=3.8
\\
39x_1+74.5x_2+24x_3&=-2.4
\\
58x_1+75.2x_2-51x_3&=12.9
\end{cases}$$

在规划求解窗口编写如下代码

[Constraint]:
-49*x1+82*x2+36*x3=14.67
   92*x1+6*x2+x3=3.8
39*x1+74.5*x2+24*x3=-2.4
58*x1+75.2*x2-51*x3=12.9

如下图,点击求解可以得到最终结果

2.2 非线性方程组求解

求解如下非线性方程组,其中\(x,y,z\)为浮点变量

$$\begin{cases}
2x^2 +3.6\sin(x)\cos(y-z)=5.8
\\
x + y^2 + 3z^2=83.58
\\
(y-x)(z-x) + x^3=10.305
\\
3xy=z+11.4
\\
x<y<z
\end{cases}$$

在规划求解窗口编写如下代码

[Constraint]:
 2*x^2 + 3.6*sin(x)*cos(y-z) = 5.8
x + y^2 + 3*z^2 = 83.58
(y-x) * (z-x) + x^3 = 10.305
3*x*y = z+11.4
x<y<z

如下图,点击求解可以得到最终结果\(x=1.50005453896144,y=3.59980165806281,z=4.80005135194315\)

2.3 线性规划求解

最大化如下目标函数,其中\(x_1,x_2,x_3,x_4\)为浮点变量

$$\begin{aligned}
最大化目标:& 0.75x_1-20x_2+0.5x_3-6x_4
\\
约束:&\begin{cases}
0.25x_1-8x_2-x_3+9x_4\leq0
\\
0.5x_1-12x_2-0.5x_3+3x_4\leq0
\\
x_3\leq1
\\
x_1\geq0,x_2\geq0,x_3\geq0,x_4\geq0
\end{cases}
\end{aligned}$$

在规划求解窗口编写如下代码

[MaxExpress]: 
0.75*x1-20*x2+0.5*x3-6*x4 

[Constraint]: 
0.25*x1-8*x2-x3+9*x4≤0 
0.5*x1-12*x2-0.5*x3+3*x4≤0 
x3≤1 
x1>=0,x2>=0,x3>=0,x4>=0

如下图,点击求解可以得到最终结果

2.4 非线性整数规划

最大化如下目标函数,其中\(x_1,x_2,x_3,x_4,x_5\)为整数变量

$$\begin{aligned}
最大化目标:&x_1^2+x_2^2+3x_3^2+4x_4^2+2x_5^2-8x_1-2x_2-3x_3-x_4-2x_5
\\
约束:&\begin{cases}
x_1+x_2+x_3+x_4+x_5\leq 400
\\
x_1+2x_2+2x_3+x_4+6x_5\leq 800
\\
2x_1+x_2+6x_3\leq 200
\\
x_3+x_4+5x_5\leq 200
\\
0\leq x_1\leq 99
\\
0\leq x_2\leq 99
\\
0\leq x_3\leq 99
\\
0\leq x_4\leq 99
\\
0\leq x_5\leq 99
\end{cases}
\end{aligned}$$

在规划求解窗口编写如下代码

[MaxExpress]: 
x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5

[IntegerVariable]:
x1,x2,x3,x4,x5

[Constraint]: 
x1+x2+x3+x4+x5<=400
x1+2*x2+2*x3+x4+6*x5<=800
2*x1+x2+6*x3<=200
x3+x4+5*x5<=200
0<=x1<=99
0<=x2<=99
0<=x3<=99
0<=x4<=99
0<=x5<=99

如下图,点击求解可以得到最终结果\(x_1=50,x_2=99,x_3=0,x_4=99,x_5=20\)

3 视频

直线与相机位姿以及观测之间的关系

约定
A. \(\mathbf{a}^T\mathbf{b}\)默认为一个数值。
B. \(\mathbf{q}_i\)表示矩阵\(Q\)的第\(i\)个行组成的向量。
C. \(q_i\)表示向量\(\mathbf{q}\)的第\(i\)个元素。
D.
$$
\mathbf{q}^\wedge =\begin{bmatrix}
0 & -q_2 & q_1\\
q_2 & 0 & -q_0\\
-q_1 & q_0 & 0
\end{bmatrix}
$$

问题
现在已知一直线,在相机平面内的一些观测点,其中相机位姿为\(R\)、\(\mathbf{t}\),相机参数为\(K\),其中直线上一个点在相机里的像素坐标为\((u,v)\),现在想建立直线和相机位姿以及观测直接的关系

这里设空间直线为如下的参数方程
(1)
$$
\mathbf{p} = \mathbf{n}k + \mathbf{m}
$$
公式中约定\(\mathbf{n}\)为单位方向向量,然后将直线上一个点往相机里投影,根据投影方程,可以得到
(2)
$$
s
\begin{bmatrix}
u\\
v\\
1
\end{bmatrix}
=
K(R\mathbf{p}+\mathbf{t})
\\
\Downarrow
\\
sK^{-1}\begin{bmatrix}
u\\
v\\
1
\end{bmatrix}
=R\mathbf{p}+\mathbf{t}
$$

这里令
(3)
$$
\mathbf{a} = K^{-1}\begin{bmatrix}
u\\
v\\
1
\end{bmatrix}
$$

将公式(3)带入(2)并约掉\(s\)有
(4)
$$
\begin{bmatrix}
a_0/a_2\\
a_1/a_2
\end{bmatrix}
=
\begin{bmatrix}
(\mathbf{r}_0 \cdot \mathbf{p} + t_0)/(\mathbf{r}_2 \cdot \mathbf{p} + t_2)\\
(\mathbf{r}_1 \cdot \mathbf{p} + t_1)/(\mathbf{r}_2 \cdot \mathbf{p} + t_2)
\end{bmatrix}
\\
\Downarrow
\\
\begin{bmatrix}
(\mathbf{r}_0 – a_0/a_2 \mathbf{r}_2)\cdot \mathbf{p}\\
(\mathbf{r}_1 – a_1/a_2 \mathbf{r}_2)\cdot \mathbf{p}
\end{bmatrix}
=
\begin{bmatrix}
t_2a_0/a_2 – t_0\\
t_2a_1/a_2 – t_1
\end{bmatrix}
$$

进一步,令
(5)
$$
V = \begin{bmatrix}
(\mathbf{r}_0 – a_0/a_2 \mathbf{r}_2)^T\\
(\mathbf{r}_1 – a_1/a_2 \mathbf{r}_2)^T
\end{bmatrix}_{2 \times 3}
\\
\mathbf{w} = \begin{bmatrix}
t_2a_0/a_2 – t_0\\
t_2a_1/a_2 – t_1
\end{bmatrix}
$$
则将公式(1),(5)带进公式(4)有
(6)
$$
\mathbf{w}=V\mathbf{p}=kV\mathbf{n}+V\mathbf{m}
$$
消去公式(6)中的\(k\)有
(7)
$$
\frac{w_0 – \mathbf{v}_0 \cdot \mathbf{m}}{w_1-\mathbf{v}_1 \cdot \mathbf{m}}=\frac{\mathbf{v}_0\cdot \mathbf{n}}{\mathbf{v}_1\cdot \mathbf{n}}
\\
\Downarrow
\\
(w_0 \mathbf{v}_1-w_1\mathbf{v}_0)\cdot \mathbf{n} +\mathbf{m}^T (\mathbf{v}_1\mathbf{v}_0^T-\mathbf{v}_0\mathbf{v}_1^T)\mathbf{n}=0
$$

为了简洁,这里记
(8)
$$
\mathbf{b}=w_0 \mathbf{v}_1-w_1\mathbf{v}_0
\\
\mathbf{c}^\wedge = \mathbf{v}_1\mathbf{v}_0^T-\mathbf{v}_0\mathbf{v}_1^T
$$
将(8)带进(7)有
(9)
$$
(\mathbf{b}^T + \mathbf{m}^T\mathbf{c}^\wedge)\mathbf{n} = \mathbf{n}^T(\mathbf{b} – \mathbf{c}^\wedge \mathbf{m}) =0
$$
由于\(\mathbf{n}\)是单位方向向量,因此这里不妨设置
(10)
$$
\mathbf{n}=\begin{bmatrix}
\sin \theta\\
\cos \theta \cos \beta\\
\cos \theta \sin \beta
\end{bmatrix}
$$
则最终带入方程有
(11)
$$
\sin \theta (b_0 + c_2m_1 – c_1m_2) +\cos \theta \cos \beta (b_1 + c_0m_2 – c_2m_0) +\cos \theta \sin \beta(b_2+c_1m_0 – c_0m_1)=0
$$
从公式(11)可以看出,已知多组相机参数以及观测,则可以拟合出最佳直线参数\(\mathbf{n}\)与\(\mathbf{m}\)。

同样地,如果已知直线与相机参数,以及观测,则要求直线上某个点的3d位置,根据公式(6)可以得到
(12)
$$
k = \frac{\mathbf{n}^TV^T(\mathbf{w}-V\mathbf{m})}{\mathbf{n}^TV^TV\mathbf{n}}
$$
最终将\(k\)代入公式(1)的直线方程,可得到3d点位置。

李代数求导

下面主要记录一些李代数涉及到的求导,其中主要提供左扰动与常规加法更新求导两种方式。因为不同的求导方式,会影响变量的更新法则,所以要注意区分。

推导过程主要用到的公式可参考李群李代数工具

1. SO(3) 求导
(1.1) 左扰动求导
$$
\begin{aligned}
\frac{ \partial R\mathbf{p}}{ \partial \phi}
&=\lim_{\delta \phi  \rightarrow 0} \frac{\exp(\delta \phi^\wedge)R\mathbf{p} – R\mathbf{p}}{\delta \phi}\\
& \approx \lim_{\delta \phi  \rightarrow 0} \frac{(I + \delta \phi^\wedge)R\mathbf{p} – R\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{\delta \phi^\wedge R\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{-(R\mathbf{p})^\wedge \delta \phi}{\delta \phi}\\
&= -(R\mathbf{p})^\wedge
\end{aligned}
$$

(1.2) 常规加法更新求导
$$
\begin{aligned}
\frac{ \partial R\mathbf{p}}{ \partial \phi}
&=\frac{ \partial \exp(\phi^\wedge)\mathbf{p}}{ \partial \phi}\\
&=\lim_{\delta \phi  \rightarrow 0} \frac{\exp((\phi +\delta \phi)^\wedge)\mathbf{p} – R\mathbf{p}}{\delta \phi}\\
& \approx \lim_{\delta \phi  \rightarrow 0} \frac{\exp((J_l \delta \phi)^\wedge)R\mathbf{p} – R\mathbf{p}}{\delta \phi}\\
& \approx \lim_{\delta \phi  \rightarrow 0} \frac{(I+(J_l \delta \phi)^\wedge)R\mathbf{p} – R\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{(J_l \delta \phi)^\wedge R\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{-(R\mathbf{p})^\wedge J_l \delta \phi}{\delta \phi}\\
&= -(R\mathbf{p})^\wedge J_l
\end{aligned}
$$

(1.3) 左扰动求导
$$
\begin{aligned}
\frac{ \partial R^T\mathbf{p}}{ \partial \phi}
&=\lim_{\delta \phi  \rightarrow 0} \frac{(\exp(\delta \phi^\wedge)R)^T\mathbf{p} – R^T\mathbf{p}}{\delta \phi}\\
&=\lim_{\delta \phi  \rightarrow 0} \frac{R^T\exp(-\delta \phi^\wedge)\mathbf{p} – R^T\mathbf{p}}{\delta \phi}\\
& \approx \lim_{\delta \phi  \rightarrow 0} \frac{R^T(I – \delta \phi^\wedge)\mathbf{p} – R^T\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{-R^T\delta \phi^\wedge \mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{R^T \mathbf{p}^\wedge \delta \phi}{\delta \phi}\\
&= R^T \mathbf{p}^\wedge
\end{aligned}
$$

(1.4) 常规加法更新求导
$$
\begin{aligned}
\frac{ \partial R^T\mathbf{p}}{ \partial \phi}
&=\frac{ \partial \exp(\phi^\wedge)^T\mathbf{p}}{ \partial \phi}\\
&=\lim_{\delta \phi  \rightarrow 0} \frac{\exp((\phi +\delta \phi)^\wedge)^T\mathbf{p} – R^T\mathbf{p}}{\delta \phi}\\
& \approx \lim_{\delta \phi  \rightarrow 0} \frac{\exp(-(J_r \delta \phi)^\wedge)R^T\mathbf{p} – R^T\mathbf{p}}{\delta \phi}\\
& \approx \lim_{\delta \phi  \rightarrow 0} \frac{(I-(J_r \delta \phi)^\wedge)R^T\mathbf{p} – R^T\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{-(J_r \delta \phi)^\wedge R^T\mathbf{p}}{\delta \phi}\\
&= \lim_{\delta \phi  \rightarrow 0} \frac{(R^T\mathbf{p})^\wedge J_r \delta \phi}{\delta \phi}\\
&= (R^T\mathbf{p})^\wedge J_r
\end{aligned}
$$
上面公式中,\(\phi\)为\(R\)对应的李代数向量,\(\mathbf{p}\)与\(R\)无关。

(1.5) 左扰动求导
$$
\begin{aligned} 
\dfrac{\partial \log(R_1R_2R_3)^\vee}{\partial \phi_2} &= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1\exp(\delta\phi_2^\wedge)R_2R_3)^\vee-\log(R_1R_2R_3)^\vee}{\delta \phi_2}
\\
 &= \lim_{\delta \phi_2 \rightarrow 0}\dfrac{\log(R_1R_2R_3\exp(((R_2R_3)^T\delta \phi_2)^\wedge))^\vee-\log(R_1R_2R_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0}\dfrac{\phi'+J_r'^{-1}(R_2R_3)^T\delta \phi_2-\phi'}{\delta \phi_2}
\\
&=J_r'^{-1}(R_2R_3)^T
\end{aligned}
$$

(1.6) 常规加法更新求导
$$
\begin{aligned} 
\dfrac{\partial \log(R_1R_2R_3)^\vee}{\partial \phi_2} &= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1\exp((\delta\phi_2+\phi_2)^\wedge)R_3)^\vee-\log(R_1R_2R_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1R_2\exp((J_{r2}\delta\phi_2)^\wedge)R_3)^\vee-\log(R_1R_2R_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1R_2R_3\exp((R_3^TJ_{r2}\delta\phi_2)^\wedge))^\vee-\log(R_1R_2R_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\phi'+J_r'^{-1}R_3^TJ_{r2}\delta\phi_2-\phi'}{\delta \phi_2}
\\
&=J_r'^{-1}R_3^TJ_{r2}
\end{aligned}
$$
上面公式(1.5)-(1.6)中,\(R_1,R_2,R_3\)相互之间无关,\(\phi'\)为\(R_1R_2R_3\)对应的李代数向量,\(J_r'\)为\(R_1R_2R_3\)对应右雅克比矩阵。

(1.7) 左扰动求导
$$
\begin{aligned} 
\dfrac{\partial \log(R_1R_2^TR_3)^\vee}{\partial \phi_2} &= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1(\exp(\delta \phi_2^\wedge)R_2)^TR_3)^\vee-\log(R_1R_2^TR_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1R_2^TR_3\exp(-(R_3^T\delta\phi_2)^\wedge)^\vee-\log(R_1R_2^TR_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\phi''-J_r''^{-1}R_3^T\delta\phi_2-\phi''}{\delta \phi_2}
\\
&=-J_r''^{-1}R_3^T
\end{aligned}
$$

(1.8) 常规加法更新求导
$$
\begin{aligned} 
\dfrac{\partial \log(R_1R_2^TR_3)^\vee}{\partial \phi_2} &= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1\exp((\delta\phi_2+\phi_2)^\wedge)^TR_3)^\vee-\log(R_1R_2^TR_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1R_2^T\exp((-J_{l2}\delta\phi_2)^\wedge)R_3)^\vee-\log(R_1R_2^TR_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\log(R_1R_2^TR_3\exp((-R_3^TJ_{l2}\delta\phi_2)^\wedge))^\vee-\log(R_1R_2^TR_3)^\vee}{\delta \phi_2}
\\
&= \lim_{\delta \phi_2 \rightarrow 0} \dfrac{\phi''-J_r''^{-1}R_3^TJ_{l2}\delta\phi_2-\phi''}{\delta \phi_2}
\\
&=-J_r''^{-1}R_3^TJ_{l2}
\end{aligned}
$$

上面公式(1.7)-(1.8)中,\(R_1,R_2,R_3\)相互之间无关,\(\phi''\)为\(R_1R_2^TR_3\)对应的李代数向量,\(J_r''\)为\(R_1R_2^TR_3\)对应右雅克比矩阵。

2. SE(3) 求导
(2.1)
$$
\mathbf{z} = R \mathbf{p} + \mathbf{t}\\
\Downarrow
\\
\begin{bmatrix}
\mathbf{z}\\
1
\end{bmatrix}
= \begin{bmatrix}
R & \mathbf{t} \\ 
\mathbf{0}^T & 1
\end{bmatrix}
\begin{bmatrix}
\mathbf{p}\\
1
\end{bmatrix}
= T\mathbf{p}'
$$

(2.2) 左扰动求导
$$
\begin{aligned}
\frac{\partial T \mathbf{p}'}{ \partial \xi}
&=\lim_{\delta \xi  \rightarrow 0} \frac{\exp(\delta \xi^\wedge)T\mathbf{p}' – T \mathbf{p}'}{\delta \xi}\\
& \approx \lim_{\delta \xi  \rightarrow 0} \frac{(I + \delta \xi^\wedge)T\mathbf{p}' – T \mathbf{p}'}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{\delta \xi^\wedge T\mathbf{p}'}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{
\begin{bmatrix}
\delta \phi^\wedge & \delta \rho \\ 
\mathbf{0}^T & 0
\end{bmatrix}
\begin{bmatrix}
R & \mathbf{t} \\ 
\mathbf{0}^T & 1
\end{bmatrix}
\begin{bmatrix}
\mathbf{p} \\ 
1
\end{bmatrix}
}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{
\begin{bmatrix}
\delta \phi^\wedge (R\mathbf{p} + \mathbf{t}) + \delta \rho \\ 
0
\end{bmatrix}
}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{
\begin{bmatrix}
-(R\mathbf{p} + \mathbf{t})^\wedge \delta \phi + \delta \rho \\ 
0
\end{bmatrix}
}{\delta \xi}\\
&=\begin{bmatrix}
I & -(R\mathbf{p} + \mathbf{t})^\wedge\\ 
0  & 0
\end{bmatrix}
\end{aligned}
\\
\Downarrow
\\
\frac{\partial (R \mathbf{p} + \mathbf{t})}{ \partial \xi}= \begin{bmatrix}
I & -(R\mathbf{p} + \mathbf{t})^\wedge
\end{bmatrix}
$$

(2.3) 常规加法更新求导
$$
\begin{aligned}
\frac{\partial T \mathbf{p}'}{ \partial \xi}
&=\lim_{\delta \xi  \rightarrow 0} \frac{\exp((\xi + \delta \xi)^\wedge)\mathbf{p}' – T \mathbf{p}'}{\delta \xi}\\
& \approx\lim_{\delta \xi  \rightarrow 0} \frac{\exp((\mathcal{J}_l \delta \xi)^\wedge)T\mathbf{p}' – T \mathbf{p}'}{\delta \xi}\\
& \approx \lim_{\delta \xi  \rightarrow 0} \frac{(I + (\mathcal{J}_l\delta \xi)^\wedge)T\mathbf{p}' – T \mathbf{p}'}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{(\mathcal{J}_l\delta \xi)^\wedge T\mathbf{p}'}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{
\begin{bmatrix}
(J_l\delta \phi)^\wedge & J_l \delta \rho + Q_l \delta \phi \\ 
\mathbf{0}^T & 0
\end{bmatrix}
\begin{bmatrix}
R & \mathbf{t} \\ 
\mathbf{0}^T & 1
\end{bmatrix}
\begin{bmatrix}
\mathbf{p} \\ 
1
\end{bmatrix}
}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{
\begin{bmatrix}
(J_l \delta \phi)^\wedge (R\mathbf{p} + \mathbf{t}) + J_l \delta \rho + Q_l \delta \phi \\ 
0
\end{bmatrix}
}{\delta \xi}\\
&= \lim_{\delta \xi  \rightarrow 0} \frac{
\begin{bmatrix}
-(R\mathbf{p} + \mathbf{t})^\wedge (J_l\delta \phi) + J_l \delta \rho + Q_l \delta \phi \\ 
0
\end{bmatrix}
}{\delta \xi}\\
&=\begin{bmatrix}
J_l & Q_l -(R\mathbf{p} + \mathbf{t})^\wedge J_l\\ 
0  & 0
\end{bmatrix}
\end{aligned}
\\
\Downarrow
\\
\frac{\partial (R \mathbf{p} + \mathbf{t})}{ \partial \xi}= \begin{bmatrix}
J_l & Q_l -(R\mathbf{p} + \mathbf{t})^\wedge J_l
\end{bmatrix}
$$

结合(1.1),(1.2),(2.2),(2.3)可以得到,对于同属于一个SE3中的变量,有如下结果

(2.4) 左扰动求导
$$
\begin{aligned}
\frac{\partial J_l \rho}{ \partial \phi} &=
\frac{\partial \mathbf{t}}{ \partial \phi} = -\mathbf{t}^\wedge
\end{aligned}
$$

(2.5) 常规加法更新求导
$$
\frac{\partial \mathbf{t}}{ \partial \phi} = Q_l -\mathbf{t}^\wedge J_l
$$

(2.4)(2.5)不仅给出了偏移与\(\phi\)的求导关系,还给出了雅可比矩阵求导与\(\phi\)的关系

(2.6) 左扰动求导
$$
\begin{aligned}
\frac{\partial \mathbf{t}}{ \partial \rho} &= I
\end{aligned}
$$

(2.7) 常规加法更新求导
$$
\frac{\partial \mathbf{t}}{ \partial \rho} = J_l
$$

从公式(2.4),(2.5)可以得到左雅克比矩阵的求导

(2.8)左扰动求导

$$
\begin{aligned}
\frac{\partial J_l \mathbf{p}}{ \partial \phi} &= -(J_l \mathbf{\mathbf{p}})^\wedge
\end{aligned}
$$

(2.9) 常规加法更新求导

$$
\begin{aligned}
\frac{\partial J_l \mathbf{p}}{ \partial \phi} &=Q_l -(J_l \mathbf{\mathbf{p}})^\wedge J_l
\end{aligned}
$$

注意,这里\(Q_l\)当中的\(\rho\)用\(\mathbf{p}\)代替。

而对于右雅克比矩阵,有

(2.10)

$$
\begin{aligned}
\frac{\partial J_r \mathbf{p}}{ \partial \phi} = \frac{\partial R^T J_l \mathbf{p}}{ \partial \phi}= \frac{\partial R^T (J_l \mathbf{p})}{ \partial \phi} + R^T \frac{\partial   J_l \mathbf{p}}{ \partial \phi}
\end{aligned}
$$

将(1.3),(1.4),(2.8),(2.9)代入公式(2.10)可以得到

(2.11)左扰动求导

$$
\frac{\partial J_r \mathbf{p}}{ \partial \phi} = \mathbf{0}
$$

(2.12)常规加法更新求导

$$
\frac{\partial J_r \mathbf{p}}{ \partial \phi} =R^T(\mathbf{p}^\wedge J_l + Q_l – (J_l\mathbf{p})^\wedge J_l)
$$

(2.13) 左扰动求导

$$
\begin{aligned}
\dfrac{\partial \log(T_1T_2T_3)^\vee}{\partial\xi_2} &\approx \lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1\exp(\delta \xi_2^\wedge)T_2T_3)^\vee – \log(T_1T_2T_3)^\vee}{\delta \xi_2}
\\
&= \lim_{\delta \xi_1 \rightarrow 0}
\dfrac{\log(T_1T_2T_3\exp((\mathcal{T}(T_2T_3)^{-1} \delta \xi_2)^\wedge))^\vee-\log(T_1T_2T_3)^\vee}{\delta \xi_2}
\\
&\approx \lim_{\delta \xi_2 \rightarrow 0}
\dfrac{\log(\exp((\xi_{123}+\mathcal{J}_{r123}^{-1}\mathcal{T}(T_2T_3)^{-1} \delta \xi_2)^\wedge))^\vee-\log(\exp(\xi_{123}^\wedge))^\vee}{\delta \xi_2}
\\
&=\mathcal{J}_{r123}^{-1}\mathcal{T}(T_2T_3)^{-1}
\end{aligned}
$$

(2.14) 常规加法求导
$$
\begin{aligned}
\dfrac{\partial \log(T_1T_2T_3)^\vee}{\partial\xi_2} &\approx \lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1\exp((\xi_2 +\delta \xi_2)^\wedge)T_3)^\vee – \log(T_1T_2T_3)^\vee}{\delta \xi_2}
\\
&\approx \lim_{\delta \xi_1 \rightarrow 0}
\dfrac{\log(T_1T_2\exp((\mathcal{J}_{r2}\delta \xi_2)^\wedge)T_3)^\vee-\log(T_1T_2T_3)^\vee}{\delta \xi_2}
\\
&= \lim_{\delta \xi_1 \rightarrow 0}
\dfrac{\log(T_1T_2T_3\exp((\mathcal{T}(T_3)^{-1}\mathcal{J}_{r2}\delta \xi_2)^\wedge))^\vee-\log(T_1T_2T_3)^\vee}{\delta \xi_2}
\\
&\approx \lim_{\delta \xi_2 \rightarrow 0}
\dfrac{\log(\exp((\xi_{123}+\mathcal{J}_{r123}^{-1}\mathcal{T}(T_3)^{-1}\mathcal{J}_{r2} \delta \xi_2)^\wedge))^\vee-\log(\exp(\xi_{123}^\wedge))^\vee}{\delta \xi_2}
\\
&=\mathcal{J}_{r123}^{-1}\mathcal{J}_{r2}\mathcal{T}(T_3)^{-1}
\end{aligned}
$$
上面公式中 \(\mathcal{J}_{r123}\)为 \(T_1T_2T_3\)的右雅可比矩阵,\(\mathcal{J}_{r2}\)为 \(T_2\)的右雅可比矩阵。

(2.15) 左扰动求导
$$
\begin{aligned}
\dfrac{\partial \log(T_1T_2^{-1}T_3)^\vee}{\partial\xi_2} &\approx \lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1(\exp(\delta \xi_2^\wedge)T_2)^{-1}T_3)^\vee – \log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&=\lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1T_2^{-1}\exp(-\delta \xi_2^\wedge)T_3)^\vee – \log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&= \lim_{\delta \xi_1 \rightarrow 0}
\dfrac{\log(T_1T_2^{-1}T_3\exp((-\mathcal{T}(T_3)^{-1} \delta \xi_2)^\wedge))^\vee-\log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&\approx \lim_{\delta \xi_2 \rightarrow 0}
\dfrac{\log(\exp((\xi_{123}'-\mathcal{J}_{r123}'^{-1}\mathcal{T}(T_3)^{-1} \delta \xi_2)^\wedge))^\vee-\log(\exp(\xi_{123}'^\wedge))^\vee}{\delta \xi_2}
\\
&=-\mathcal{J}_{r123}'^{-1}\mathcal{T}(T_3)^{-1}
\end{aligned}
$$

(2.16) 常规加法求导
$$
\begin{aligned}
\dfrac{\partial \log(T_1T_2^{-1}T_3)^\vee}{\partial\xi_2} &\approx \lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1(\exp((\xi_2+\delta \xi_2)^\wedge)T_2)^{-1}T_3)^\vee – \log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&\approx \lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1(T_2\exp((\mathcal{J}_{r2}\delta \xi_2)^\wedge))^{-1}T_3)^\vee – \log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&= \lim_{\delta \xi_2 \rightarrow 0} \dfrac{\log(T_1\exp(-(\mathcal{J}_{r2}\delta \xi_2)^\wedge))T_2^{-1}T_3)^\vee – \log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&= \lim_{\delta \xi_1 \rightarrow 0}
\dfrac{\log(T_1T_2^{-1}T_3\exp((-\mathcal{T}(T_2^{-1}T_3)^{-1}\mathcal{J}_{r2} \delta \xi_2)^\wedge))^\vee-\log(T_1T_2^{-1}T_3)^\vee}{\delta \xi_2}
\\
&\approx \lim_{\delta \xi_2 \rightarrow 0}
\dfrac{\log(\exp((\xi_{123}'-\mathcal{J}_{r123}'^{-1}\mathcal{T}(T_2^{-1}T_3)^{-1}\mathcal{J}_{r2} \delta \xi_2)^\wedge))^\vee-\log(\exp(\xi_{123}'^\wedge))^\vee}{\delta \xi_2}
\\
&=-\mathcal{J}_{r123}'^{-1}\mathcal{T}(T_2^{-1}T_3)^{-1}\mathcal{J}_{r2}
\end{aligned}
$$
上面公式中 \(\mathcal{J}_{r123}'\)为 \(T_1T_2^{-1}T_3\)的右雅可比矩阵,\(\mathcal{J}_{r2}\)为 \(T_2\)的右雅可比矩阵

 

(2.17) 左扰动求导

$$\begin{aligned}
\dfrac{ \partial T^{-1}\mathbf{p}'}{ \partial \xi}&=\lim_{\delta \xi \rightarrow 0}\dfrac{(\exp(\delta \xi^\wedge)T)^{-1}\mathbf{p}'-T^{-1}\mathbf{p}'}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0}\dfrac{T^{-1}\exp(-\delta \xi^\wedge)\mathbf{p}'-T^{-1}\mathbf{p}'}{\delta \xi}
\\
&\approx \lim_{\delta \xi \rightarrow 0}\dfrac{T^{-1}(I-\delta \xi^\wedge)\mathbf{p}'-T^{-1}\mathbf{p}'}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0}\dfrac{-T^{-1}\delta \xi^\wedge\mathbf{p}'}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0}-\dfrac{
\begin{bmatrix}
R^T&-R^T\mathbf{t}
\\
\mathbf{0}^T&1
\end{bmatrix}
\begin{bmatrix}
\delta \phi^\wedge & \delta \rho
\\
\mathbf{0}^T & 0
\end{bmatrix}
\begin{bmatrix}
\mathbf{p}
\\
1
\end{bmatrix}
}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0}-\dfrac{
\begin{bmatrix}
R^T(\delta \phi^\wedge \mathbf{p} + \delta \rho)
\\
0
\end{bmatrix}
}{\delta \xi}
\\
&=\begin{bmatrix}
-R^T & R^T\mathbf{p}^\wedge
\\
0&0
\end{bmatrix}
\end{aligned}$$

(2.18)左扰动求导

$$\begin{aligned}
\dfrac{\partial ((T^{-1})^T\mathbf{p}')}{\partial \xi} &= \lim_{\delta \xi \rightarrow 0} \dfrac{((\exp(\delta \xi^\wedge)T)^{-1})^T\mathbf{p}'-(T^{-1})^T\mathbf{p}'}{\delta \xi}
\\
&=
\lim_{\delta \xi \rightarrow 0} \dfrac{\exp(-\delta \xi^\wedge)^T(T^{-1})^T\mathbf{p}'-(T^{-1})^T\mathbf{p}'}{\delta \xi}
\\
&\approx \lim_{\delta \xi \rightarrow 0} \dfrac{(I-\delta \xi^\wedge)^T(T^{-1})^T\mathbf{p}'-(T^{-1})^T\mathbf{p}'}{\delta \xi}
\\
&= \lim_{\delta \xi \rightarrow 0} \dfrac{(-\delta \xi^\wedge)^T(T^{-1})^T\mathbf{p}'}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0} \dfrac{
\begin{bmatrix}
(-\delta \phi^\wedge)^T &\mathbf{0}
\\
-\delta \rho^T & 0
\end{bmatrix}
\begin{bmatrix}
R & \mathbf{0}
\\
-\mathbf{t}^TR & 1
\end{bmatrix}
\begin{bmatrix}
\mathbf{p}
\\
1
\end{bmatrix}
}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0}\dfrac{\begin{bmatrix}
\delta \phi^\wedge R\mathbf{p}
\\
-\delta \rho^T R\mathbf{p}
\end{bmatrix}
}{\delta \xi}
\\
&=\lim_{\delta \xi \rightarrow 0}\dfrac{\begin{bmatrix}
– (R\mathbf{p})^\wedge \delta \phi
\\
– (R\mathbf{p})^T\delta \rho
\end{bmatrix}
}{\delta \xi}
\\
&=\begin{bmatrix}
\mathbf{I} \times 0 & -(R \mathbf{p})^\wedge
\\
-(R\mathbf{p})^T & \mathbf{0}^T
\end{bmatrix}
\end{aligned}$$

(2.19)左扰动求导

$$\begin{aligned} 
\dfrac{\partial (T^T\mathbf{p}')}{\partial \xi} &=\lim_{\delta \xi\rightarrow 0}\dfrac{(\exp(\delta \xi^\wedge)T)^T\mathbf{p}'-T^T\mathbf{p}'}{\delta \xi}
\\
&=\lim_{\delta \xi\rightarrow 0}\dfrac{T^T(
    I+\delta \xi^\wedge)^T\mathbf{p}'-T^T\mathbf{p}'}{\delta \xi}
\\
&=\lim_{\delta \xi\rightarrow 0}\dfrac{
\begin{bmatrix}
R^T & \mathbf{0}
\\
\mathbf{t}^T & 1
\end{bmatrix}
\begin{bmatrix}
-\delta \phi^\wedge & \mathbf{0}
\\
\delta \rho^T & 0
\end{bmatrix}
\begin{bmatrix}
\mathbf{p}
\\
1
\end{bmatrix}
}{\delta \xi}
\\
&=\lim_{\delta \xi\rightarrow 0}\dfrac{
\begin{bmatrix}
R^T\mathbf{p}^\wedge \delta\phi
\\
\mathbf{t}^T\mathbf{p}^\wedge\delta \phi+\mathbf{p}^T\delta \rho
\end{bmatrix}
}{\delta \xi}
\\
&=\begin{bmatrix}
\mathbf{I} \times 0 & R^T\mathbf{p}^\wedge
\\
\mathbf{p}^T & \mathbf{t}^T\mathbf{p}^\wedge
\end{bmatrix}
\end{aligned}$$

线性EKF

PS:EKF本身是一种数学工具,具体的原理推导可以参考文献[1]。对于使用者来说,只要理解这个工具对应的数学模型,以及对应系数的更新方式,就可以套用EKF的求解模式进行求解。

1 模型
在离散线性系统中,对于某种状态\(\mathcal{x}\)

定义运动方程 
(1.1)
$$
\mathcal{x}_{k} = A_{k-1} \mathcal{x}_{k-1} + \mathcal{v}_k + \mathcal{w}_k
$$

公式中\(\mathcal{w}_k \in \mathcal{N}(\mathbf{0}, Q_k)\)的过程噪声,其中\(A_{k-1}\)称为转移矩阵,\(\mathcal{x}_0 \in \mathcal{N}(\mathbf{0}, P_0)\),\(P_0\)为\(\mathcal{x}_0\)的初始协方差矩阵。

同样,定义观测方程
(1.2)
$$
\mathcal{y}_{k} = C_k \mathcal{x}_k + \mathcal{n}_k
$$
公式中\(\mathcal{n}_k \in \mathcal{N}(\mathbf{0}, R_k)\)的测量噪声,其中\(C_{k}\)称为观测矩阵。

 

2 卡尔曼滤波
下面为常规的扩展卡尔曼滤波更新过程[1]
1.预测 
(2.1)
$$
\begin{aligned} 
P_k &= A_{k-1}P_{k-1}A_{k-1}^T + Q_k\\
\mathcal{x}_k &=A_{k-1}\mathcal{x}_{k-1}+\mathcal{v}_k
\end{aligned} 
$$

2.卡尔曼增益
(2.2)
$$
\begin{aligned}
K &= P_kC_k^T(C_kP_kC_k^T+R_k)^{-1}\\
\end{aligned}
$$

3.更新
(2.3)
$$
\begin{aligned}
P_k&= (I – KC_k)P_k\\
\mathcal{x}_k &= \mathcal{x}_k + K(\mathcal{y}_k – C_k\mathcal{x}_k)
\end{aligned}
$$

即,每次可以使用公式(2.1),(2.2),(2.3)进行过程更新。

3 理解
在每一轮更新中,先根据运动方程进行预测,然后结合观测方程估计出一个修正系数,最后使用系数对变量进行修正。

仔细观察上面更新过程以及对应数学模型,可以发现,在本质上来说可以通过运动方程直接递推出下一个状态,但是考虑到一些噪音(为了提高精度),这里会使用观测方程来对递推出的状态进行修正,也因此可以理解成两类方程耦合

 

注意,上面公式中协方差矩阵\(P_{k}\)本应该服务的对象,下面将以SLAM当中的例子加以说明。

4 SLAM当中的EKF 

4.1 运动方程
对于位姿\(T\),以及它对应的广义速度\(\varpi\),它们满足运动学方程  
(4.1)
$$
\dot{T} = \varpi ^\wedge T
$$
对上面的两个变量,分边添加一个微小的扰动,然后忽略二阶小量有 
(4.2)
$$
\begin{aligned} 
\frac{d}{dt}((I + \delta \xi^\wedge)T) &= (\delta \varpi + \varpi)^\wedge ((I +\delta \xi ^\wedge) T)\\
\delta \dot{\xi}^\wedge T + \dot{T} + \delta \xi ^\wedge \dot{T} &=\varpi ^\wedge T + \varpi ^\wedge \delta \xi^\wedge T + \delta \varpi^\wedge T + \delta \varpi^\wedge \delta \xi ^\wedge T\\
& \approx \varpi ^\wedge T + \varpi ^\wedge \delta \xi^\wedge T + \delta \varpi^\wedge T
\end{aligned}
$$

将(4.1)带入(4.2)有 
(4.3)
$$
\begin{aligned}
\delta \dot{\xi}^\wedge &= \varpi ^\wedge \delta \xi ^\wedge – \delta \xi ^\wedge \varpi ^\wedge + \delta \varpi ^\wedge\\
&=(\varpi ^\curlywedge \delta \xi )^\wedge+ \delta \varpi ^\wedge\\
&=(\varpi ^\curlywedge \delta \xi + \delta \varpi)^\wedge
\end{aligned}
\\
\Downarrow
\\
\delta \dot{\xi} = \varpi ^\curlywedge \delta \xi + \delta \varpi
$$
一般把(4.1)称为标称运动学方程,(4.2)称为扰动运动学方程。
对上面公式进行离散化可得到 
(4.4)
$$
\begin{aligned}
T_k &= \exp((t_k – t_{k-1}) \varpi ^\wedge)T_{k-1}\\
\delta \xi_k &= \exp((t_k – t_{k-1})\varpi ^\curlywedge) \delta \xi _{k-1} + \mathcal{w}_k
\end{aligned}
$$

\(t_{k}\)表示对应时刻值,\(\mathcal{w}_{k} \in \mathcal{N}(0, Q_k)\)为过程噪声。 

4.2 观测方程
这里假设具有如下的观测模型 
(4.5)
$$
\mathcal{y} =  T \mathcal{p} + n
$$
同样,两边同时加一个微小的扰动有 
(4.6)
$$
\begin{aligned}
\mathcal{y} + \delta \mathcal{y} &= (I+ \delta \xi^\wedge)Tp + n\\
&= Tp + \delta \xi ^\wedge Tp + n\\
&= Tp + (Tp)^ {\odot} \delta \xi + n
\end{aligned}
$$
这里将(4.5)的标称解带入(4.6)有 
(4.7)
$$
\delta y = (Tp)^{\odot} \delta \xi + n
$$
即有最终离散化公式 
(4.8)
$$
\delta y_k = (Tp)^{\odot} \delta \xi_k + n_k
$$

上面\(\mathcal{n}_k \in \mathcal{N}(\mathbf{0}, R_k)\)


4.3 EKF过程
这里先约定如下变量 
(4.9)
$$
\begin{aligned}
O_k &= \exp((t_k – t_{k-1}) \varpi_{k-1} ^\wedge) \\
G_k &= \exp((t_k – t_{k-1})\varpi_{k-1} ^\curlywedge)\\
C_k &= (T_kp_k)^{\odot}
\end{aligned}
$$

最终结果如下 
4.3.1 预测   
(4.10)
$$
\begin{aligned}
T_k &= O_k T_{k-1}\\
P_k &= G_k P_{k-1}G_k^T + Q_k
\end{aligned}
$$

4.3.2 卡尔曼增益 
(4.11)
$$
\begin{aligned}
K = P_k C_k^T (C_k P_k C_k^T + R_k)^{-1}
\end{aligned}
$$

4.3.3 矫正 
(4.12)
$$
\begin{aligned}
P_k &= (I – KC_k)P)k\\
\mathcal{y}_k' &= T_k \mathcal{p}_k\\
T_k &= \exp((K(\mathcal{y}_k – \mathcal{y}_k'))^\wedge)T_k
\end{aligned}
$$
 参考
[1] 高翔,谢晓佳. 机器人学中的状态估计[M]. 西安交通大学出版社,2018

位姿的刚体变换与坐标系转换

1 前言
在一些场景里,可能会遇到坐标系转换与刚体变换的问题,对于空间点来说,它们的变换形式是一致的,而对于带旋转的位姿来说,它们的变换方式却不一致。

一个对象进行刚体变换,整体表现为旋转与平移,不改变坐标系本身约定的性质(各个轴之间的相对关系不会改变)。刚体变换可以理解成在同一个约定的坐标系下,从一个相对坐标系到另一个相对坐标系的变换过程!(PS:世界坐标系也可看成一个特殊的相对坐标系)

坐标系转换,从名字上就可以看出,其本身是对约定的系统进行改变。比如,左右坐标系的转换,它们各个轴的相对关系因坐标系的不同而不同。

不管是刚体变换还是坐标系之间的转换,主要的难点是旋转矩阵的表现形式。本文这里仅简单记录。

注意: 
(1) 本文中没有特别说明的,其表示方式均为世界坐标系的表示方式。 
(2) 本文当中的3D点\(\mathbf{p}\)为了表示方便,均以齐次坐标形式表示,即 
(1.1)
$$
\mathbf{p}=\begin{bmatrix} 
x\\
y\\
z\\
1
\end{bmatrix}
$$

 

2 刚体变换

刚体变换的主体思想是,对操作对象进行旋转与平移,这里约定有如下(2.1)的变换矩阵,然后用这个矩阵对操作对象进行操作即达到刚体变换效果。 
(2.1)
$$
T =  \begin{bmatrix} 
R & \mathbf{t}\\
\mathbf{0}^T & 1
\end{bmatrix}
$$

 

2.1 点的变换
对于空间点\(\mathbf{p}_0\),假设在刚体变换\(T\)后变为\(\mathbf{p}_1\),那么变换前后的关系可以表示为 
(2.2)

$$
\mathbf{p}_1 = T \mathbf{p}_0
$$

 

2.2 位姿的变换
对于空间中,有一个相机的位姿为\(T_0\),假设在刚体变换\(T\)后变为\(T_1\)。同样,空间点\(\mathbf{p}_0\),在刚体变换\(T\)后变为\(\mathbf{p}_1\)。由于刚体变换不会改变空间中对象的相对位置信息,因此可知点\(\mathbf{p}_0\)在相机\(T_0\)里的相对位置不会改变,因此可以得到 
(2.3)

$$
T_0\mathbf{p}_0 = T_1\mathbf{p}_1=T_1T\mathbf{p}_0
$$

由于对任意\(\mathbf{p}_0\)均要满足公式(2.3),因此有 
(2.4)

$$
T_1 = T_0 T^{-1}
$$

注意,这里的位姿\(T_0、T_1\)是世界坐标系到相机坐标系的表示方式, 即如果一个世界坐标系的点\(\mathbf{p}_0\),通过公式\(T_0 \mathbf{p}_0\)将变换到相机坐标系下的位置。如果\(T_0、T_1\)相机坐标系到世界坐标系的表示方式, 即如果一个世界坐标系的点\(\mathbf{p}_0\),通过公式\(T_0^{-1} \mathbf{p}_0\)将变换到相机坐标系下的位置,那么根据公式(2.4),可以得到如下位姿变换结果

(2.5)

$$
T_1 = TT_0
$$

3 坐标系变换
有时我们会在左右手坐标系之间进行转换,有时我们需要将x、y轴进行对调,有时我们需要将当前坐标系按某种形式进行旋转……这一系列的坐标系变换,可以看成一个镜像、旋转、偏移的组合,这里以如下的变换矩阵来表示 
(3.1)
$$
T =  \begin{bmatrix} 
R' & \mathbf{t}\\
\mathbf{0}^T & 1
\end{bmatrix}
$$
注意:这里的\(R'\)的行列式不一定为1,也可能为-1。

 

3.1 点的坐标系变换
对于空间点\(\mathbf{p}_0\),假设在坐标系变换\(T\)后变为\(\mathbf{p}_1\),那么变换前后的关系可以表示为 
(3.2)

$$
\mathbf{p}_1 = T \mathbf{p}_0
$$

 

3.2 位姿的坐标系变换

3.2.1 相机坐标系按约定变换,世界坐标系按约定变换
对于空间中,有一个相机的位姿为\(T_0\),假设经过坐标系变换系数\(T\)后变为\(T_1\)。同样,空间点\(\mathbf{p}_0\),在经过坐标系变换\(T\)后变为\(\mathbf{p}_1\)。由于\(\mathbf{p}_0\)在相机\(T_0\)里的相对位置不会改变。但是,需要注意的是,虽然相对位置不会改变,但是按约定,坐标系进行了变化(例如,在世界坐标系,我们约定x轴向前,y轴向右,z轴向上;在相机坐标系,我们也会默认x轴向前,y轴向右,z轴向上),则可得到如下结果 

这里假设,在变换前的坐标系下,\(\mathbf{p}_0\)相对于\(T_0\)的位置为\(\mathbf{p}'\),则可以得到 
(3.3)

$$
\mathbf{p}' = T_0 \mathbf{p}_0
$$

\(\mathbf{p}'\)在经过坐标系变换\(T\)后变为\(\mathbf{p}''\),则有  
(3.4)

$$
\mathbf{p}'' = T \mathbf{p}'
$$

由于在变换后的坐标系下,\(\mathbf{p}_1\)相对于\(T_1\)的位置也为\(\mathbf{p}''\),即 
(3.5)

$$
\mathbf{p}'' = T_1 \mathbf{p}_1
$$
根据公式(3.2)-(3.5)可得到 
(3.6)

$$
(T_1T – TT_0)\mathbf{p}_0 = \mathbf{0}
$$
由于对任意\(\mathbf{p}_0\)均要满足公式(3.6),因此有 
(3.7)

$$
T_1 = T T_0 T^{-1}
$$

对(3.7)进行展开有   
(3.8)

$$
\begin{bmatrix} 
R_1 & \mathbf{t}_1\\
\mathbf{0}^T & 1
\end{bmatrix}
=
\begin{bmatrix} 
RR_0R^{-1} &R\mathbf{t}_0 + \mathbf{t} – RR_0R^{-1}\mathbf{t}\\
\mathbf{0}^T & 1
\end{bmatrix}
$$

即旋转部分与偏移部分的变换结果分别为 
(3.9)

$$
\begin{aligned} 
R_1 &= RR_0R^{-1}\\
\mathbf{t}_1 &= R\mathbf{t}_0 + \mathbf{t} – RR_0R^{-1}\mathbf{t}
\end{aligned} 
$$

特别地,当\(\mathbf{t}=\mathbf{0}\)时,有 
(3.10)
$$
\begin{aligned} 
R_1 &= RR_0R^{-1}\\
\mathbf{t}_1 &= R\mathbf{t}_0
\end{aligned} 
$$

 

3.2.2 相机坐标系不变,世界坐标系按约定变换

这里同样使用3.2.1当中的符号表达,由于这里约定相机坐标系不变,因此(3.5)式变为

(3.11)

$$
\mathbf{p}' = T_1 \mathbf{p}_1
$$

由公式(3.2),(3.3),(3.11)可以得到

(3.12)

$$
(T_1T – T_0)\mathbf{p}_0 = \mathbf{0}
$$

由于对任意\(\mathbf{p}_0\)均要满足公式(3.12),因此有 
(3.13)

$$
T_1 = T_0T^{-1}
$$

 

3.2.3 相机坐标系按约定变换,世界坐标系不变

这里同样使用3.2.1当中的符号表达,由于这里约定世界坐标系不变,因此(3.5)式变为

(3.14)

$$
\mathbf{p}'' = T_1 \mathbf{p}_0
$$

由公式(3.3),(3.4),(3.14)可以得到

(3.15)

$$
(TT_0 – T_1)\mathbf{p}_0 = \mathbf{0}
$$

由于对任意\(\mathbf{p}_0\)均要满足公式(3.15),因此有 
(3.16)

$$
T_1 = TT_0
$$

 

3.3 相机坐标系变换总结

最终,相机坐标系变换方式,可分3种情况进行:

情况1相机坐标系按约定变换,世界坐标系按约定变换,这时使用公式(3.7)

情况2相机坐标系不变,世界坐标系按约定变换,这时使用公式(3.13)

情况3相机坐标系按约定变换,世界坐标系不变,这时使用公式(3.16)

 

4 例子
现在有一个原始的相机位姿为\(T_0\) 
(4.1)
$$
T_0 =  \begin{bmatrix} 
R_0 & \mathbf{t}_0\\
\mathbf{0}^T & 1
\end{bmatrix}
$$
因为一些接口需要,我们希望将原始坐标系的x轴指向新坐标系的z轴反向,原坐标系的y轴指向新坐标系的x轴反向,原坐标系的z轴指向新坐标系y轴。根据需求,可得到变换矩阵的旋转部分 
(4.2)

$$
R = \begin{bmatrix} 
0 & -1 & 0\\
0 & 0 & 1\\
-1 & 0 & 0
\end{bmatrix}
$$
根据需求,偏移部分为 
(4.3)

$$
\mathbf{t} = \begin{bmatrix} 
0\\
0\\
0
\end{bmatrix}
$$
因此最终得到扩展后的变换矩阵 
(4.4)
$$
T = \begin{bmatrix} 
0 & -1 & 0 & 0\\
0 & 0 & 1 &0 \\
-1 & 0 & 0 &0\\
0 & 0 & 0 & 1
\end{bmatrix}
$$

按公式(3.7),可得到最终的坐标系变换结果 
$$
T_1 = T T_0 T^{-1}
$$

 

5 延伸

5.1 按指定向量旋转

很多时候,我们可能希望得到某些点以\(\mathbf{p}_0\)为原心,由向量\(\mathbf{n}_1\) 旋转到向量\(\mathbf{n}_2\)后的位置,这个本质上是求公式(3.2)当中的变换矩阵\(T\)。

很明显,第一步,可以将所有点以\(\mathbf{p}_0\)为原心,进行平移;第二步,由于向量\(\mathbf{n}_1\) 与向量\(\mathbf{n}_2\)叉乘后的单位轴向量\(\mathbf{a}\)同时垂直于\(\mathbf{n}_1\)、\(\mathbf{n}_2\),因此\(\mathbf{a}\)可以作为旋转轴,而旋转角度,可以通过余弦公式计算得出。因此可先按如下方式确定对应变量:

(5.1)

$$
\begin{aligned}
\mathbf{a} &=\frac{ \mathbf{n}_1 \times \mathbf{n}_2}{|\mathbf{n}_1||\mathbf{n}_2|}\\
\theta &= \mathbf{arccos}(\frac{\mathbf{n}_1 \cdot \mathbf{n}_2}{|\mathbf{n}_1||\mathbf{n}_2|})
\end{aligned}
$$

结合SLAM之李群李代数工具中的公式(1.5),可以得到旋转矩阵\(R\),最终变换矩阵形式如下

(5.2)

$$
T = \begin{bmatrix}  
R & -R\mathbf{p}_0\\
\mathbf{0} & 1
\end{bmatrix}
$$

SO(3)与Sim(3)插值

在常规的插值中,我们会根据n个已知点的值,插值出给定点的结果。其中主要利用的思想是,通过这n个点建立某种数学模型,进而推导出这个待插值点的结果。本文这里仅考虑2个已知点的情况,对应到我们常规插值中,就表现为线性插值。但本文这里,由于插值出的结果不是一个数值,是一个需要满足某种性质的矩阵,因此,不能单纯采用矩阵的线性加减。

1 SO(3)旋转矩阵插值

这里直接给出一种插值模型

(1)$$\begin{aligned}  
R &= (R_1R_0^T)^aR_0\\   
&= R_{10}^aR_0\\   
&=\exp(\varphi ^ \wedge)^aR_0\\ 
&=\exp(a \varphi ^ \wedge)R_0 \\
&=\exp((a \varphi )^ \wedge)R_0 
\end{aligned}$$

上面公式中的\(a\)表示插值点的值,其取值范围在[0,1],当\(a = 0\)时,对应值为旋转矩阵\(R_0\),当\(a = 1\)时,对应值为旋转矩阵\(R_1\),公式中的\(\varphi\)对应的矩阵为 \(R_1R_0^T\)

 

2 Sim(3)位姿矩阵插值

先申明如下SE(3)矩阵格式

(2)$$
T_0 = \begin{bmatrix}  
s_0R_0 & t_0\\
0 & 1
\end{bmatrix}\\T_1 = \begin{bmatrix}  
s_1R_1 & t_1\\
0 & 1
\end{bmatrix}\\T = \begin{bmatrix}  
sR & t\\
0 & 1
\end{bmatrix}$$

即,现在我们想找到一种插值模型,使得\(a = 0\)时,插值的\(T = T_0\);当\(a = 1\)时,\(T = T_1\),当\(a\)取其它值时,插值的\(T\)满足Sim(3)性质。其中插值点\(a\)取值在[0,1]。

即现在需要通过\(T_0,T_1\)插值出\(T\)当中的\(s,R,t\)

2.1 无光心约束插值

与旋转矩阵插值类似,可以构建如下的插值模型

(3)$$\begin{aligned}  
T &= (T_1T_0^{-1})^aT_0\\   
&= T_{10}^aT_0\\   
&=\exp(\xi ^ \wedge)^aT_0\\ 
&=\exp(a \xi ^ \wedge)T_0 \\
&=\exp((a \xi) ^ \wedge)T_0
\end{aligned}$$

这里将公式(3)展开有

(4)$$T=\begin{bmatrix}  
(s_1/s_0)^a(R_1R_0^T)^a & a(t_1 – (s_1/s_0)R_1R_0^Tt_0)\\
0 & 1
\end{bmatrix}
\begin{bmatrix}  
s_0R_0 & t_0\\
0 & 1
\end{bmatrix}$$

通过公式(4)可以得到最终的插值结果

(5)$$\begin{aligned}  
s &= s_0(\frac{s_1}{s_0})^a\\
R &= (R_1R_0^T)^aR_0\\
t &= (\frac{s_1}{s_0})^a(R_1R_0^T)^at_0 + a(t_1-\frac{s_1}{s_0}R_1R_0^Tt_0)
\end{aligned}$$

2.2 有光心约束插值

当我们使用Sim(3)来表示一个位姿的时候,我们知道一种位姿在空间会对应一个相机中心坐标。当给定两个位姿,然后在它们之间进行插值时,我们可以假设它们插值的相机中心也满足两个位姿相机中心的线性插值。即有如下约束

(6)$$\frac{R^Tt}{s} =(1-a) \frac{R_0^Tt_0}{s_0}  + a\frac{R_1^Tt_1}{s_1}$$

公式(6)是一个含有3个等式的约束,如果插值结果的旋转\(R\)与缩放部分\(s\)采用公式(5)的结果,那结合公式(6)刚好可以求出插值的\(t\),因此最终的插值结果可表示为

(7)$$\begin{aligned}  
s &= s_0(\frac{s_1}{s_0})^a\\
R &= (R_1R_0^T)^aR_0\\
t &=sR ((1-a) \frac{R_0^Tt_0}{s_0}  + a\frac{R_1^Tt_1}{s_1})
\end{aligned}$$

复数的基本运算

PS:这篇文章来自于2011年我写在网易博客上的内容,由于网易博客现已关闭,因此特把此文章搬迁到这里,供有缘人参考,指正!

提醒:当你某一天发现同一个复数领域内的算式表达得出不同的结果的时候,请不要太过怀疑,因为这就是复数的多值性!

已知:

1

这里设\(x\)、\(y\)、\(m\)、\(n\)为实数,\(i\)表示虚数单位,即\(i = \sqrt{-1}\)

$$
\begin{aligned}  
A &= x + yi\\   
&= e^{a + \psi i}\\   
&=|A|e^{\psi i}\\ 
&=|A|(\cos \psi + \sin \psi i)
\end{aligned}
$$

$$
\begin{aligned}  
B &= m + ni\\   
&= e^{b + \beta i}\\   
&=|B|e^{\beta i}\\ 
&=|B|(\cos \beta + sin \beta i)
\end{aligned}
$$

其中

$$
|A| = e^a = \sqrt{x^2+y^2}, a = \ln(|A|), \psi =\arctan(y/x)\\
|B| = e^b = \sqrt{m^2+n^2}, b = \ln(|B|), \beta = \arctan(n/m)
$$

以\(A\)为例,\(|A|\)叫做复数\(A\)的模,\(\psi\)为复数\(A\)的复角,一般把\(-\pi < \psi < \pi\)的值叫做主值,此时定义\(arg(A)=\psi\),当\(A \neq 0\)的时候,主值由如下方式确定

(1)、当\(x > 0\)时,\(\psi = \arctan(y/x)\)

(2)、当\(x=0,y>0\)时,\(\psi =\frac{\pi}{2}\)

(3)、当\(x=0,y<0\)时,\(\psi =-\frac{\pi}{2}\)

(4)、当\(x<0,y \geqslant 0\)时,\(\psi =\arctan(y/x)+\pi\)

(5)、当\(x<0,y<0\)时,\(\psi =\arctan(y/x)-\pi\)

基本运算:

注意,下面的最后的结果都会变成\(x +yi\)的形式)

$$
A+B=(x+m)+(y+n)i\\
A-B=(x-m)+(y-n)i
$$


$$
\begin{aligned}  
AB &= (x+yi)(m+ni)\\
&= xm+xni+ymi+ynii\\ 
&=(xm-yn)+(xn+ym)i
\end{aligned}
$$


$$
\begin{aligned}  
\frac{A}{B} &= \frac{x+yi}{m+ni}\\
&= \frac{(x+yi)(m-ni)}{(m+mi)(m-ni)}\\
&=\frac{(xm+yn)+(ym-xn)}{m^2+n^2}\\
&=\frac{xm+yn}{m^2+n^2}+\frac{ym-xn}{m^2+n^2}i
\end{aligned}
$$

$$
\begin{aligned}  
\log(A) &=\log(|A|e^{\psi i})\\
&=\log(|A|)+\log(e^{\psi i})\\
&=\frac{\log(x^2+y^2)}{2}+\psi i\\
&=\log(|A|)+\mathbf{arg}(A)i
\end{aligned}
$$


$$
\begin{aligned}  
A^B &= (e^{a+\psi i})^{m+ni}\\
&=e^{(a+\psi i)(m+ni)}\\
&=e^{(am-\psi n)+(an+\psi m)i}\\
&=e^{(am-\psi n)}e^{(an+\psi m)i}\\
&=e^{(am-\psi n)}(\cos(an+\psi m)+\sin(an+\psi m)i)\\
&=e^{(am-\psi n)}\cos(an+\psi m)+e^{(am-\psi n)}\sin(an+\psi m)i
\end{aligned}
$$

对于\(A^{0.5} = \sqrt(A)\)这种特殊情况,由于开方在一些计算中调用频繁,因此不建议直接使用上面求次方的方式进行计算,而建议直接使用下面提到的算法。

设\(B=\sqrt{A}\),即\(A=B^2\)可以通过对于表达式最终反求出结果,这里就不贴过程,直接给出答案

$$
\sqrt{A}=\sqrt{\frac{\sqrt{x^2+y^2}+x}{2}} + \mathbf{sign}(y)\sqrt{\frac{\sqrt{x^2+y^2}-x}{2}}i
$$

注意,本文在在上面公式中,当\(y<0\)时,\(\mathbf{sign}(y)=-1\);当\(y=0\)时,\(\mathbf{sign}(y)=0\);当\(y>0\)时,\(\mathbf{sign}(y)=1\)

$$
\sin(A)=\frac{e^{Ai}-e^{-Ai}}{2i}
$$

$$
\cos(A)=\frac{e^{Ai}+e^{-Ai}}{2}
$$

$$
\tan(A)= \frac{\sin(A)}{\cos(A)}
$$

 

$$
\cot(A)= \frac{\cos(A)}{\sin(A)}
$$

$$
\sec(A)=\frac{1}{\cos(A)}
$$

$$
\csc(A)=\frac{1}{\sin(A)}
$$

$$
\sinh(A)=\frac{e^A-e^{-A}}{2}
$$

$$
\cosh(A)=\frac{e^A+e^{-A}}{2}
$$

$$
\tanh(A)=\frac{\sinh(A)}{\cosh(A)}
$$

$$
\coth(A)=\frac{\cosh(A)}{\sinh(A)}
$$

$$
\ln_B(A)=\frac{\ln(A)}{\ln(B)}
$$

$$
\begin{aligned}  
\arcsin(A)&=-\log(Ai+\sqrt{1-A^2})i
\\
\\
\arccos(A)&=-\log(A+\sqrt{A^2-1})i
\\
\\
\arctan(A)&=-\dfrac{1}{2}\log(\dfrac{1+Ai}{1-Ai})i
\\
\\
\mathbf{arccot}(A)&=-\dfrac{1}{2}\log(\dfrac{Ai-1}{Ai+1})i
\\
\\
\mathbf{arcsinh}(A)&=\log(A+\sqrt{1+A^2})
\\
\\
\mathbf{arccosh}(A)&=\log(A+\sqrt{A^2-1})
\\
\\
\mathbf{arctanh}(A)&=\dfrac{1}{2}\log(\dfrac{1+A}{1-A})
\\
\\
\mathbf{arccoth}(A)&=\dfrac{1}{2}\log(\dfrac{A+1}{A-1})
\end{aligned}
$$

地震来了,公式更新延后

离散数据插值预测

一 前言

数据插值在很多地方都会有所应用。很多时候,我们需要根据已知点(姑且称为点吧)去插值(预测)出我们想要点对应的数据值。有时一个点可能只有1个变量(一维插值问题),有时可能有多个变量(多维插值问题)。这时如何去比较真实的预测出对应的值,这就是问题所在。

 

二 求解

1、提到插值,我们首先想到的是数值分析书上介绍的一些内容,几乎你在每一本数值分析书上都能找到关于插值的介绍。上面主要针对的是一维插值,因为书里介绍得比较多,这里不再详叙!

2、这里我们需要解决的是多维问题的插值,而我们更多的了解却是基于一维的插值算法。因为一维问题比较简单,那么我们是否可以像多维非线性拟合等问题求解将多维问题转化为一维问题来进行求解呢。答案当然是肯定的。

3、一般地,我们认为两个点比较靠近时,它们所赋予的物理量(需要插值出的数据变量)的值比较接近,因为向量范数具有距离的概念且它可将多维变量转换成一维变量,因此我们可利用范数来降低求解维度。为了编写程序方便,可将多维变量根据范数定义类似采用如下算式
(1)$$r=(x_1^2+x_2^2+……)^v$$
算式中\(x_1,x_2,……\)为多维自变量,\(v\)为一种控制距离的参数,\(r\)即为转换成的一维自变量。

4、权重思想插值

假如我们现在已知点的数据为\(m\)个,我们如何根据这\(m\)个数据插值出指定点的数据呢。其中一个最简单的思想是权重思想。

想象一下,已知一个正方形四个顶点对应的值,那如果想求形心处的值,一般我们会将4个顶点值相加求和,最后将和除以4,其实对每个点来说,它的权重系数为\(\frac{1}{4}\),求形心点的值时是每个点的值乘以对应权重系数,最后求和所得。

而如何构造每个点对应的权重系数则成了一个关键点。因为第3点已经求得每个点距离插值点的距离\(r\),而我们认为两个点比较靠近时,它们所赋予的物理量(需要插值出的数据变量)的值比较接近,即相近的点影响较大。这时很自然可按照如下算式来定义其权重
(2)$$ K(r)=\frac{1}{r}$$

可以使用公式(2)来构造一个权重系数产生函数。当然我们可以认为距离\(r\)越近影响越大,且在靠近一定程度影响比较大,远离一段距离后会很快衰减,则可考虑如下类似的算式(3)(注意,你也可以按其它方式进行构造)
(3)$$K(r)=\frac{exp(-ar)}{R}$$

\(a\)为控制系数,\(r\)是相对距离,\(R\)是实际距离。
则最后每个点的权重系数可按如下算式(4)求得
(4)$$k_i=\frac{K(r_i)}{∑K(r_j)},(j=1,2,……,m)$$

如果每个点对应的值为\(f\),则最后可得到预测点的\(f\)值为算式(5)
(5)$$f=k_1f_1+k_2f_2+……+k_mf_m$$

5、径向基函数插值
这种插值算法插值出的效果比较厉害,个人理解其就好比在一个\(n\)(插值的自变量个数为\(n\))维空间里织了一个很大的网,然后在这个网里来插值出最佳值。简单理解,它就是假设函数值为算式(6)
(6)
$$f(x)=c_1K(r_1)+c_2K(r_2)+……+c_1K(r_m)$$

算式中\(r_i\)表示\(x-x_i\)通过公式(1)获得,\(c_i\)表示求解的系数。

因为它有\(m\)个未知数\(c\)需要求解,但是在这个函数里,它必须满足每个点的值,即可将已知的\(m\)个点带进去进而求解得到每个系数\(c\),最后将求解点\(x\)带进去求解即可。

对于(6)算式当中的函数\(K\)一般有一定要求,具体可参看相关文献,这里对于紧支型的函数\(K\),对于\(r \leqslant 1\)有
文献[1] 

(7)$$
\begin{align*}
K(r)&=1-r 
\\
K(r)&=(1-r)^3(3r+1) 
\\
K(r)&=(1-r)^5(8r^2+5r+1) 
\\
K(r)&=(1-r)^2 
\\
K(r)&=(1-r)^4(4r+1) 
\\
K(r)&=(1-r)^6(35r^2+18r+3) 
\\
K(r)&=(1-r)^8(32r^3+25r^2+8r+1) 
\\
K(r)&=(1-r)^3 
\\
K(r)&=(1-r)^5(5r+1) 
\\
K(r)&=(1-r)^7(16r^2+7r+1)
\end{align*} 
$$
文献[2] 

(8)$$
\begin{align*}
K(r)&=(1-r)^2(2+r) 
\\
K(r)&=(1-r)^3(1+3r+r^2) 
\\
K(r)&=(1-r)^3(8+9r+3r^2) 
\\
K(r)&=(1-r)^4(4+16r+12r^2+3r^3) 
\\
K(r)&=(1-r)^5(1+5r+9r^2+5r^3+r^4) 
\\
K(r)&=(1-r)^4(16+29r+20r^2+5r^3) 
\\
K(r)&=(1-r)^5(8+40r+48r^2+25r^3+5r^4) 
\\
K(r)&=(1-r)^6(6+36r+82r^2+72r^3+30r^4+5r^5) 
\\
K(r)&=(1-r)^7(5+35r+101r^2+147r^3+101r^4+35r^5+5r^6)
\end{align*}
$$

6、其它细节问题

6.1、插值点的确定

考虑到一些实际问题,很多时候我们往往只需要进行局部插值即可。即我们需要确定哪些点作为已知点来进行插值。文献当中往往以插值点为圆心,一个半径\(R\)来确定出作为插值的已知点。如果数据分布很散乱且在我们不知道如何确定\(R\)的情况下(即我们不知道数据空间分布规律),选出一个合适的\(R\)有一定困难。这时建议可以一个数据量\(N\)来划定局部区域范围,即每次插值时只选择距离插值点最近的\(N\)个点进行插值即可。

6.2、半径\(r\)规范化

由于不同的数据处理,按(1)式计算出的每个距离大小无法确定。为了便于统一,可将所有距离\(r\)进行规范化。如算式(9)。这时的\(r\)可直接带进公式(2)、(3)以及公式(7)、(8)的紧支函数\(K\)里。

(9)$$r_i=\frac{r_i}{r_0}$$
\(r_0\)为所有半径当中的最大值.

当然,我在编写径向紧支函数求解问题时,我并不是直接采用公式(9),而是采用公式(10)
(10)$$r_i=\frac{r_i}{br_0}$$

\(r_0\)为动态距离(不一定是最大值),\(b\)为一个控制参数一般取16,公式当中的\(r_i\)均不会大于1,否则重新选取一个合适的\(r_0\)。

6.3、数据再插值的后处理

对于一些插值结果,可以进行简单的在加工。比如,在进行二维正方形网格插值时,现在已经插值出了各个网格节点值,则可再根据网格节点值插值出各个正方形网格型心对应的值,然后再根据各个网格型心的值再对各个节点进行简单插值,得到平滑效果。

参考文献
[1] Wendland H.Piecewise polynomial,positive definite and compactly supported radial functions of minimal degree[J]. Advances in computational Mathematics,1995,4(1):389-396.  

[2] Wu Z. Compactly supported positive definite radial functions[J]. Advances in Computational Mathematics, 1995, 4(1): 283-292.

位姿与平面耦合

1 前言
在使用单目相机的无人驾驶中,有时我们需要构建位姿与现实中某一个平面的关系,以便为后续其它操作提供基础!

 

2 原理

2.1 位姿与平面间的耦合

假设有相机1(\(s_1\)、\(R_1\)、\(t_1\))与相机2(\(s_2\)、\(R_2\)、\(t_2\)),一个平面在相机1坐标系下的方程为
(1)
$$
n^Tp+d=0
$$
上面公式中\(n\)是相机平面法向量,\(p\)是一个在相机1坐标系下路平面上的点,\(d\)实际就是相机到平面的距离。对公式(1)进行变形有
(2)
$$
1=-\frac{n^Tp}{d}
$$

根据SLAM之位姿到位姿的相对测量中的公式(5)可以得到相机2与相机1的相对位姿关系
(3)
$$
\begin{aligned}
R &= R_2R_1^Ts_2/s_1
\\
t &= t_2 – Rt_1
\\
s_2R_2 &= s_1RR_1
\\
t_2 &= t + Rt_1
\end{aligned}
$$

现在(1)的平面上,有一个点被相机1与相机2看到,其在世界坐标系下的位置为\(p_w\),其在相机1坐标系下的坐标为\(p_c\),其在相机1平面上看到的坐标为\(c_1(u_1,v_1,1)\),在相机2平面上看到点坐标为\(c_2(u_2,v_2,1)\),假设相机1与相机2的内参分别为\(K_1\)、\(K_2\),现在令
(4)
$$
p_1=K_1^{-1}c_1\\
p_2=K_2^{-1}c_2
$$

根据针孔相机模型可以得到
(5)
$$
z_1c_1 = K_1(s_1R_1p_w+t_1)=K_1p_c\\
\Downarrow \\
p_c = z_1K_1^{-1}c_1 = z_1p_1
$$

(6)
$$
\begin{aligned}
z_2c_2&=K_2(s_2R_2p_w+t_2)
\\
&=K_2(s_1RR_1p_w+t+Rt_1)
\\
&=K_2(R(s_1R_1p_w+t1)+t)
\\
&=K_2(Rp_c+t)
\\
&=K_2(Rp_c+t(-\frac{n^Tp_c}{d}))
\\
&=K_2(R-\frac{tn^T}{d})p_c
\\
&=K_2(R-\frac{tn^T}{d})z_1p_1
\end{aligned}
\\
\Downarrow
\\
z_2p_2=z_2K_2^{-1}c_2=(R-\frac{tn^T}{d})z_1p_1
$$

公式(6)建立了平面与位姿的关系。这里需要理解:

1)公式当中是2个位姿共同耦合1个平面。

2)耦合的方程中还包含两个相机共同观测到平面上的点。

3)公式当中的平面是相机1坐标系下的平面,如果要建立世界坐标系下平面与位姿的关系,则将平面方程进行坐标系转换一下即可(相机1坐标系转到世界坐标系)

 

2.2 根据观测求取相对测量与平面
公式(5)、(6)当中的\(z_1\)、\(z_2\)是使得相机归一化后的最后一维的值,这里令\(s = z_2/z_1\)则公式(6)最终可表示为 
(7)
$$
sp_2=(R-\frac{tn^T}{d})p_1
$$

且复原到\(c_1\)、\(c_2\)有 
(8)
$$
sc_2=K_2(R-\frac{tn^T}{d})K_1^{-1}c_1=Hc_1
$$
SLAM当中将\(H\)叫做单应矩阵,一般会使用这个矩阵来恢复出相对测量\(R\)、\(t\),下面将公式(8)展开有
(9)
$$
s\begin{bmatrix}
u_2\\ 
v_2\\
1
\end{bmatrix}
=
\begin{bmatrix}
h_1 & h_2 & h_3\\
h_4 & h_5 & h_6\\
h_7 & h_8 & h_9
\end{bmatrix}
\begin{bmatrix}
u_1\\ 
v_1\\
1
\end{bmatrix}
\\
\Downarrow
\\
\begin{aligned}
su_2 &= h_1u_1+h_2v_1+h_3\\
sv_2 &= h_4u_1+h_5v_1+h_6\\
s&=h_7u_1+h_8v_1+h_9
\end{aligned}
$$

在公式(9)当中消去\(s\)有
(10)
$$
\begin{aligned}
u_2 &= \frac{h_1u_1+h_2v_1+h_3}{h_7u_1+h_8v_1+h_9}\\
v_2 &= \frac{h_4u_1+h_5v_1+h_6}{h_7u_1+h_8v_1+h_9}
\end{aligned}
$$
仔细观察公式(10)可以发现,对于\(H\)中的元素同时乘以一个非0的常数,并不改变公式(10)的结果,且这样操作后等价于\(H\)的特征值同时乘以这个非0常数。因此这里不妨假设\(h_9=1\),当然也可以假设\(H\)的其它元素值为1(这里需要注意,假设的元素的值本身可能为0,所以保险做法是取几个元素置为1分别尝试计算\(R\)、\(t\)),因此最终方程(10)可变为
(11)
$$
\begin{bmatrix}
u_1 & v_1 & 1 & 0 & 0 & 0 & -u_1u_2 & -v_1u_2
\\
0 & 0 & 0 & u_1 & v_1 & 1 & -u_1v_2 & -v_1v_2
\end{bmatrix}
\begin{bmatrix}
h_1\\
h_2\\
h_3\\
h_4\\
h_5\\
h_6\\
h_7\\
h_8
\end{bmatrix}
=
\begin{bmatrix}
u_2\\
v_2
\end{bmatrix}
$$

当有多组观测点都落在这个平面时,可以通过组装的线性方程组求解出一个最佳的\(H\),由于在求解过程中假设\(h_9=1\),因此真实的\(H\)可表示为\(\alpha H\)

因此根据公式(8)最终可以得到
(12)
$$
(R-\frac{tn^T}{d}) = \alpha K_2^{-1}HK_1 = \alpha G
$$

对于公式(1)的平面方程来说,两边同时乘以一个非0系数,不改变平面性质,因此对于公式(12)当中的\(d\),不妨令其等于1(实际当中一般需要建立相机与平面的关系时,这个相机一般距离平面有一定距离,因此参数\(d\)为0的可能性很小),则公式(12)则变成
(13)
$$
(R-tn^T) = \alpha G
$$
在上面的公式中\(G\)已经得到,对于公式(13)的形式,根据Ezio Malis[1]的结果,最终可以求解出\(R\)、\(t\)、\(n\)、\(\alpha \)

 

2.3 已知相对测量求平面方程

和上面类似,仔细观察公式(7)可知,同时缩放\(n\)、\(d\)一个非0系数,它们的平面性质不会改变(公式(1)当中两边同时乘以一个非0常数,表示的平面不变),这时不妨设$d=1$,则公式(7)最终变为
(14)
$$
sp_2=Rp_1-tn^Tp_1=Rp_1-tp_1^Tn
$$
将公式(14)中的未知数都移动到一边有
(15)
$$
tp_1^Tn+sp_2=Rp_1
$$

从公式(15)可以看出,这个里面有3个线性方程4个未知数(\(n\)为3个未知数,\(s\)为1个未知数)。

对于一个平面上的多个匹配点对,由于每个(匹配对)点对应的\(s\)不一样,因此如果平面上同时有\(m\)个点被两个相机同时看见,这时根据(15)可以构建一个$3m \times (3 + m)$的线性矩阵来求解出平面方程系数\(n\)。

如果两个相机的连线与所在平面平行,现在设两个相机连线在世界坐标系下的方向向量为\(v\),则它们与公式(1)当中的\(n\)满足如下关系
(16)
$$
n^TR_1v=0
$$
上面公式(16)是一个3个未知数(\(n\))的线性方程,因此如果将公式(16)带入公式(15),\(n\)的求解维度可以进一步降低到2维。进而可以进行求解。

参考

[1] Ezio Malis, Manuel Vargas. Deeper understanding of the homography decomposition for vision-based control. 2007