采用逆深度参数表达的BA问题导数推导
由于大部分的slam算法均采用逆深度参数表达地图点的结构,但是网上对该方法的介绍比较少因此本文将详细说明其推导过程。逆深度参数表达具有优化变量少、能表达非常远的点以及分布接近高斯分布等优势,这也是大家选择逆深度参数的原因。
图结构表达
从上图可以看出一个约束项Ep12E_{p_{12}}Ep12连接三个节点分别是当前帧、参考帧和逆深度参数表达的地图点。这是因为在逆深度参数的表达下需要当前帧、参考帧位姿及地图点才能共同确定误差函数,这将在下文中看到。
优化节点
优化节点分为两类,关键帧的相机位姿和地图点坐标。其中关键帧的相机位姿维护SE(3)表达的状态,而地图点则是逆深度参数。
相机位姿
状态量采用四元数+平移向量的方式表达,即[q|t],对于微小的角度更新方法如下:
tk+1=tk+δtqk+1=qk⊕δq
t_{k+1} = t_{k} + \delta_t \\
q_{k+1} = q_{k} \oplus \delta_{q}
tk+1=tk+δtqk+1=qk⊕δq
其中⊕\oplus⊕代表的是四元数更新,一般采用四阶Runge-Kunta法更新四元数。
∂pf∂δξ=(TP)⊙=[I−pf∧0T0T]
\frac{\partial p_f}{\partial \delta_{\xi}} = (TP)^{\odot}=\left[\begin{matrix}
I & -p_f^{\wedge} \\
0^T & 0^T
\end{matrix}\right]
∂δξ∂pf=(TP)⊙=[I0T−pf∧0T]
其中pfp_fpf为当前坐标系下的XYZXYZXYZ表达地图点,导数的具体求解过程参考附录。
逆深度参数地图点
逆深度参数采用MSCKF表示:
pf=[xfyfzf]=1ρ[αβ1]
p_f = \left[\begin{matrix}
x_f \\
y_f \\
z_f
\end{matrix}\right] = \frac{1}{\rho}\left[\begin{matrix}
\alpha \\
\beta \\
1
\end{matrix}\right]
pf=⎣⎡xfyfzf⎦⎤=ρ1⎣⎡αβ1⎦⎤
λ=[αβρ]=[xf/zfyf/zf1/zf]
\lambda = \left[\begin{matrix}
\alpha \\
\beta \\
\rho
\end{matrix}\right]= \left[\begin{matrix}
x_f/z_f \\
y_f/z_f \\
1/z_f
\end{matrix}\right]
λ=⎣⎡αβρ⎦⎤=⎣⎡xf/zfyf/zf1/zf⎦⎤
∂pf∂λ=[1ρ0−αρ201ρ−βρ200−1ρ2]
\frac{\partial{p_f}}{\partial \lambda} = \left[\begin{matrix}
\frac{1}{\rho} & 0 & -\frac{\alpha}{\rho^2} \\
0 & \frac{1}{\rho} & -\frac{\beta}{\rho^2}\\
0 & 0 & -\frac{1}{\rho^2}
\end{matrix}\right]
∂λ∂pf=⎣⎢⎡ρ1000ρ10−ρ2α−ρ2β−ρ21⎦⎥⎤
其中pfp_fpf为参考坐标系下地图点的三维坐标,xf,yf,zfx_f,y_f,z_fxf,yf,zf分别为对应数值,λ\lambdaλ为逆深度参数表达的坐标而二者的导数关系也一并给出。推导公式参考附录。
边约束项
误差计算的思路在于把参考帧下的地图点转换到当前帧坐标系下,在投影到当前帧像素平面上与对应的特征点坐标求差值。因此误差函数为:
e=[umvm]−[uv]
e = \left[\begin{matrix}
u_m\\v_m
\end{matrix}\right]−\left[\begin{matrix}
u\\v
\end{matrix}\right]
e=[umvm]−[uv]
[umvm]=H(Twc−1TwfG(λ))
\left[\begin{matrix}
u_m\\v_m
\end{matrix}\right]= H(T_{wc}^{-1}T_{wf}G(\lambda))
[umvm]=H(Twc−1TwfG(λ))
其中[umvm]\left[\begin{matrix}u_m\\v_m\end{matrix}\right][umvm]是重投影的点,HHH函数根据相机模型将相机坐标系下的点投影到像素坐标系中,而GGG函数则将逆深度参数转换为参考帧下的3d坐标,对于Twc,TwfT_{wc},T_{wf}Twc,Twf则分别是世界系下当前相机帧的位姿和参考帧的相机位姿。
更进一步,将误差函数展开成为更加具体的形式。
[umvm]=[fxczc+cxfyczc+cy] \left[\begin{matrix} u_m \\ v_m \end{matrix}\right] = \left[\begin{matrix} f \frac{x_c}{z_c} + c_x \\ f \frac{y_c}{z_c} + c_y \end{matrix}\right] [umvm]=[fzcxc+cxfzcyc+cy]
pc=[xcyczc]=Rwc−1pw−Rwc−1twc p_c=\left[\begin{matrix} x_c \\ y_c \\ z_c \end{matrix}\right] = R_{wc}^{-1}p_w - R_{wc}^{-1}t_{wc} pc=⎣⎡xcyczc⎦⎤=Rwc−1pw−Rwc−1twc
pw=Rwfpf+twf p_w = R_{wf}p_f + t_{wf} pw=Rwfpf+twf
pf=1ρ[αβ1] p_f = \frac{1}{\rho}\left[\begin{matrix} \alpha \\ \beta \\ 1 \end{matrix}\right] pf=ρ1⎣⎡αβ1⎦⎤
将误差函数通过一步步转换成为逆深度参数坐标,在求导时便可轻易的采用链式法则进行计算。
对地图点的导数
∂e∂λ=∂e∂pc∂pc∂pf∂pf∂λ \frac{\partial e}{ \partial \lambda} = \frac{\partial e}{ \partial p_c} \frac{\partial p_c}{ \partial p_f} \frac{\partial p_f}{ \partial \lambda} ∂λ∂e=∂pc∂e∂pf∂pc∂λ∂pf
由于∂pf∂λ\frac{\partial p_f}{ \partial \lambda}∂λ∂pf已知,因此只需要求解∂e∂pc,∂pc∂pf\frac{\partial e}{ \partial p_c} ,\frac{\partial p_c}{ \partial p_f}∂pc∂e,∂pf∂pc
∂e∂pc=−[∂u∂xc∂u∂yc∂u∂zc∂v∂xc∂v∂yc∂v∂zc]=−[fzc0−fxczc20fzc−fyczc2] \frac{\partial e}{ \partial p_c} = -\left[\begin{matrix} \frac{\partial u}{ \partial x_c} & \frac{\partial u}{ \partial y_c} & \frac{\partial u}{ \partial z_c} \\ \frac{\partial v}{ \partial x_c} & \frac{\partial v}{ \partial y_c} & \frac{\partial v}{ \partial z_c} \end{matrix}\right] = - \left[\begin{matrix} \frac{f}{z_c} & 0 & -\frac{f x_c}{z_c^2} \\ 0 & \frac{f}{ z_c} & -\frac{f y_c}{z_c^2} \end{matrix}\right] ∂pc∂e=−[∂xc∂u∂xc∂v∂yc∂u∂yc∂v∂zc∂u∂zc∂v]=−[zcf00zcf−zc2fxc−zc2fyc]
∂pc∂pf=Rcf \frac{\partial p_c}{ \partial p_f} = R_{cf} ∂pf∂pc=Rcf
对当前帧的导数
∂e∂δξc=∂e∂pc∂pc∂δξc
\frac{\partial e}{ \partial \delta_{\xi_c}} = \frac{\partial e}{ \partial p_c} \frac{\partial p_c}{\partial \delta_{\xi_c}}
∂δξc∂e=∂pc∂e∂δξc∂pc
其中,∂e∂pc\frac{\partial e}{ \partial p_c}∂pc∂e在对点求导时已求得,而∂pc∂δξc\frac{\partial p_c}{\partial \delta_{\xi_c}}∂δξc∂pc只要带入可得
∂pc∂δξc=[I−pc∧0T0T]
\frac{\partial p_c}{\partial \delta_{\xi_c}} = \left[\begin{matrix}
I & -p_c^{\wedge} \\
0^T & 0^T
\end{matrix}\right]
∂δξc∂pc=[I0T−pc∧0T]
在计算时通常会直接删除下面为零的部分,以便和前面矩阵维度一致。
对参考帧的导数
∂e∂δξf=∂e∂pf∂pf∂δξf
\frac{\partial e}{ \partial \delta_{\xi_f}} = \frac{\partial e}{ \partial p_f} \frac{\partial p_f}{\partial \delta_{\xi_f}}
∂δξf∂e=∂pf∂e∂δξf∂pf
同样,∂e∂pf\frac{\partial e}{ \partial p_f}∂pf∂e在对点求导时已经求得了
∂pc∂δξf=[I−pf∧0T0T]
\frac{\partial p_c}{\partial \delta_{\xi_f}} =
\left[\begin{matrix}
I & -p_f^{\wedge} \\
0^T & 0^T
\end{matrix}\right]
∂δξf∂pc=[I0T−pf∧0T]
附录
向量的求导
∂[ab]∂[xy]=[∂a∂x∂a∂y∂b∂x∂b∂y] \frac{\partial\left[\begin{matrix} a \\ b \end{matrix}\right]}{\partial \left[\begin{matrix} x \\ y \end{matrix}\right]} = \left[\begin{matrix} \frac{\partial a}{\partial x} & \frac{\partial a}{\partial y}\\ \frac{\partial b}{\partial x} & \frac{\partial b}{\partial y} \end{matrix}\right] ∂[xy]∂[ab]=[∂x∂a∂x∂b∂y∂a∂y∂b]