OpenVINS 官方文档 第二部分

参考链接:OpenVINS https://docs.openvins.com/index.html


4.Measurement Update Derivations

4.1 最小均方误差 (MMSE) 估计
  考虑以下静态状态估计问题:给定 n n n 维高斯随机向量 x ∽ N ( x ^ ㊀ , P x x ㊀ ) x\backsim \mathcal{N}(\hat x^{㊀},P^{㊀}_{xx}) xN(x^,Pxx) 的先验分布(概率密度函数或pdf),和一个新的 m m m 维具有状态无关的零均值白高斯噪声 n ∽ N ( 0 , R ) n\backsim \mathcal{N}(0,R) nN(0,R) 的观测 z m = z + n = h ( x ) + n z_m = z+n=h(x)+n zm=z+n=h(x)+n,我们想要计算后验PDF的前两个(中心)时刻 p ( x ∣ z m ) p(x|z_m) p(xzm)。通常(给定非线性测量模型),我们将后验pdf近似为: p ( x ∣ z m ) ≈ N ( x ^ ⊕ , P x x ⊕ ) p(x|z_m) \approx \mathcal{N}(\hat x^{\oplus},P^{\oplus}_{xx}) p(xzm)N(x^,Pxx)。根据设计,这是MMSE估计问题的(近似)解决方案 Kay 1993

4.2 条件概率分布
  为此,我们采用贝叶斯规则:
p ( x ∣ z m ) = p ( x , z m ) p ( z m ) p(x|z_m)={p(x,z_m) \over p(z_m)} p(xzm)=p(zm)p(x,zm)
  通常,如果不施加简化的假设,则无法分析计算此条件pdf。对于手头的问题,我们首先近似 p ( z m ) ≈ N ( z ^ , P z z ) p(z_m)\approx \mathcal{N}(\hat{z},P_{zz}) p(zm)N(z^,Pzz)(如果确实如此),然后有以下联合高斯pdf(注意高斯pdf的联合仍是符合高斯的):
p ( x , z m ) = N ( [ x ^ ㊀ z ] , [ P x x P x z P z x P z z ] ) = : N ( y ^ , P y y ) p(x,z_m)=\mathcal{N}\begin{pmatrix} \left[ \begin{array}{ccc} \hat x^{㊀}\\z \\ \end{array} \right],\left[ \begin{array}{ccc} P_{xx} & P_{xz} \\ P_{zx} & P_{zz}\\ \end{array} \right] \\ \end{pmatrix}=: \mathcal{N}(\hat y, P_{yy}) p(x,zm)=N([x^z],[PxxPzxPxzPzz])=:N(y^,Pyy)
将这两个高斯代入第一个方程得到以下条件高斯pdf:

p ( x ∣ z m ) ≈ N ( y ^ , P y y ) N ( z ^ , P z z ) p(x|z_m) \approx{ \mathcal{N}(\hat y,P_{yy} ) \over \mathcal{N}(\hat z,P_{zz} )} p(xzm)N(z^,Pzz)N(y^,Pyy) = 1 ( 2 π ) n ∣ P y y ∣ / ∣ P z z ∣ e − 1 2 [ ( y − y ^ ) ⊤ P y y − 1 ( y − y ^ ) − ( z m − z ^ ) ⊤ P z z − 1 ( z m − z ^ ) ] ={1 \over \sqrt{(2 \pi)^n|P_{yy}|/|P_{zz}|}}e^{-{1 \over 2}[(y - \hat y)^{\top}P_{yy}^{-1}(y-\hat y)-(z_m-\hat z)^\top P^{-1}_{zz}(z_m-\hat z)]} =(2π)nPyy∣/∣Pzz 1e21[(yy^)Pyy1(yy^)(zmz^)Pzz1(zmz^)] = : N ( x ^ ⊕ , P x x ⊕ ) =:\mathcal{N}(\hat x^{\oplus},P^{\oplus}_{xx}) =:N(x^,Pxx)
  我们现在推导出条件均值,协方差可以计算如下:首先,我们简化分母项 ∣ P y y ∣ / ∣ P z z ∣ |P_{yy}|/|P_{zz}| Pyy∣/∣Pzz 以找到条件协方差。
∣ P y y ∣ = ∣ [ P x x P x z P z x P z z ] ∣ = ∣ P x x − P x z P z z − 1 P z x ∣ ∣ P z z ∣ |P_{yy} |= \begin{vmatrix}\left[ \begin{array}{ccc} P_{xx} & P_{xz} \\ P_{zx} & P_{zz}\\ \end{array} \right] \\ \end{vmatrix} = \begin{vmatrix} P_{xx} -P_{xz}P^{-1}_{zz}P_{zx} \\ \end{vmatrix} \begin{vmatrix} P_{zz} \\ \end{vmatrix} Pyy= [PxxPzxPxzPzz] = PxxPxzPzz1Pzx Pzz
  我们假设 P z z P_{zz} Pzz 是可逆的,并采用了舒尔补的行列式性质。因此,我们有:
∣ P y y ∣ ∣ P z z ∣ = ∣ P x x − P x z P z z − 1 P z x ∣ ∣ P z z ∣ = ∣ P x x − P x z P z z − 1 P z x ∣ {|P_{yy}| \over |P_{zz}|} = {\begin{vmatrix} P_{xx} -P_{xz}P^{-1}_{zz}P_{zx} \\ \end{vmatrix} \over |P_{zz}|} = \begin{vmatrix} P_{xx} -P_{xz}P^{-1}_{zz}P_{zx} \\ \end{vmatrix} PzzPyy=Pzz PxxPxzPzz1Pzx = PxxPxzPzz1Pzx

接下来,通过定义误差状态 r x = x − x ^ ㊀ r_x = x - \hat x^{㊀} rx=xx^ r z = z m − z ^ r_z = z_m -\hat z rz=zmz^ r y = y − y ^ r_y = y - \hat y ry=yy^ 并使用矩阵内引理,我们重写指数项如下:
( y − y ^ ) ⊤ P y y − 1 ( y − y ^ ) − ( z , − z ^ ) ⊤ P z z − 1 ( z m − z ^ ) (y - \hat y)^{\top}P_{yy}^{-1}(y-\hat y)-(z_,-\hat z)^\top P^{-1}_{zz}(z_m-\hat z) (yy^)Pyy1(yy^)(z,z^)Pzz1(zmz^) = r y ⊤ P y y − 1 r y − r z ⊤ P z z − 1 r z =r_y^{\top}P_{yy}^{-1}r_y-r_z^\top P^{-1}_{zz}r_z =ryPyy1ryrzPzz1rz = [ r x r z ] ⊤ [ P x x P x z P z x P z z ] − 1 [ r x r z ] − r z ⊤ P z z − 1 r z =\left[ \begin{array}{ccc} r_{x} \\ r{z}\\ \end{array} \right]^\top \left[ \begin{array}{ccc} P_{xx} & P_{xz} \\ P_{zx} & P_{zz}\\ \end{array} \right] ^{-1} \left[ \begin{array}{ccc} r_{x} \\ r{z}\\ \end{array} \right] - r_z^\top P^{-1}_{zz}r_z =[rxrz][PxxPzxPxzPzz]1[rxrz]rzPzz1rz = ( r x − P x z P z z − 1 r z ) ⊤ ( P x x − P x z P z z − 1 P z x ) − 1 ( r x − P x z P z z − 1 r z ) =(r_x - P_{xz}P_{zz}^{-1}r_z)^\top(P_{xx} - P_{xz}P^{-1}_{zz}P_{zx})^{-1}(r_x - P_{xz}P^{-1}_{zz}r_z) =(rxPxzPzz1rz)(PxxPxzPzz1Pzx)1(rxPxzPzz1rz)
  其中, ( p z z − 1 ) ⊤ = P z z − 1 (p_{zz}^{-1} )^\top = P_{zz}^{-1} (pzz1)=Pzz1 由于协方差矩阵是对称的。到目前为止,我们现在可以构建条件高斯pdf,如下所示:
p ( x k ∣ z m ) = 1 ( 2 π ) n ∣ P x x − P x z P z z − 1 P z x ∣ × e x p ( − 1 2 [ ( r x − P x z P z z − 1 r z ) ⊤ ( P x x − P x z P z z − 1 P z x ) − 1 ( r x − P x z P z z − 1 r z ) ] ) p(x_k|z_m)= {1 \over \sqrt{(2\pi)^n|P_{xx}-P_{xz}P_{zz}^{-1}P_{zx}|}}\times exp \begin{pmatrix} -{1\over 2} [(r_x - P_{xz}P_{zz}^{-1}r_z)^\top(P_{xx} - P_{xz}P^{-1}_{zz}P_{zx})^{-1}(r_x - P_{xz}P^{-1}_{zz}r_z)]\\ \end{pmatrix} p(xkzm)=(2π)nPxxPxzPzz1Pzx 1×exp(21[(rxPxzPzz1rz)(PxxPxzPzz1Pzx)1(rxPxzPzz1rz)])

这引出我们正在寻找的以下条件均值和协方差:
x ^ ⊕ = x ^ ㊀ + P x z P z z − 1 ( z m − z ^ ) \hat x^{\oplus} = \hat x^{㊀}+P_{xz}P_{zz}^{-1}(z_{m}-\hat z) x^=x^+PxzPzz1(zmz^) P x x ⊕ = P x x ㊀ − P x z P z z − 1 P z x P_{xx}^{\oplus} = P_{xx}^{㊀}-P_{xz}P_{zz}^{-1}P_{zx} Pxx=PxxPxzPzz1Pzx
这些是(线性)状态估计的基本方程。

4.3 线性测量更新
  作为一个特例,我们考虑一个简单的线性测量模型来说明线性MMSE估计器:
z m , k = H k x k + n k z_{m,k} = H_kx_k + n_k zm,k=Hkxk+nk z ^ k : = E [ z m , k ] = E [ H k x k + n k ] = H k x ^ k ㊀ \hat z_k:=\mathbb{E}[z_{m,k}] = \mathbb{E}[H_kx_k+n_k] = H_k\hat x_k^{㊀} z^k:=E[zm,k]=E[Hkxk+nk]=Hkx^k
有了这个,我们可以推导出协方差和相关矩阵如下:
P z z = E [ ( z m , k − z ^ k ) ( z m , k − z ^ k ) ⊤ ] P_{zz} = \mathbb{E}[(z_{m,k}-\hat z_k)(z_{m,k}-\hat z_k)^{\top}] Pzz=E[(zm,kz^k)(zm,kz^k)] = E [ ( H k x k + n k − H k x ^ k ㊀ ) ( H k x k + n k − H k x ^ k ㊀ ) ⊤ ] = \mathbb{E}[(H_kx_k+n_k-H_k\hat x_k^{㊀})(H_kx_k + n_k-H_k \hat x_k^{㊀})^\top] =E[(Hkxk+nkHkx^k)(Hkxk+nkHkx^k)] = H k P x x ㊀ H k ⊤ + R k =H_kP_{xx}^{㊀}H_{k}^\top+R_k =HkPxxHk+Rk
其中 R k R_k Rk 是离散测量噪声矩阵, H k H_k Hk 是将状态映射到测量域的测量雅可比矩阵, P x x ㊀ P_{xx}^{㊀} Pxx是当前状态协方差。
P x z = E [ ( x k − x ^ k ㊀ ) ( z m , k − z ^ k ) ⊤ ] P_{xz} = \mathbb{E}[(x_{k}-\hat x^{㊀}_k)(z_{m,k}-\hat z_k)^{\top}] Pxz=E[(xkx^k)(zm,kz^k)] = E [ ( x k − x ^ k ㊀ ) ( H k x k + n k − H k x ^ k ㊀ ) ⊤ ] = \mathbb{E}[(x_{k}-\hat x^{㊀}_k)(H_kx_k + n_k-H_k \hat x_k^{㊀})^\top] =E[(xkx^k)(Hkxk+nkHkx^k)] = P x x ㊀ H k ⊤ =P_{xx}^{㊀}H_{k}^\top =PxxHk
我们采用了噪声独立于状态的事实。将这些量代入基本方程可得到以下更新方程:
x ^ k ⊕ = x ^ k ㊀ + P k ㊀ H k ⊤ ( H k P k ㊀ H k ⊤ + R k ) − 1 ( z m , k − z ^ k ) \hat x_k^{\oplus} = \hat x^{㊀}_{k} + P_{k}^{㊀}H_k^\top(H_kP_k^{㊀}H_k^\top+R_k)^{-1}(z_{m,k}-\hat z_k) x^k=x^k+PkHk(HkPkHk+Rk)1(zm,kz^k) = x ^ k ㊀ + K r z =\hat x_k^{㊀} +Kr_z =x^k+Krz
P x x ⊕ = P k ㊀ + P k ㊀ H k ⊤ ( H k P k ㊀ H k ⊤ + R k ) − 1 ( P K ㊀ H k ⊤ ) ⊤ P_{xx}^{\oplus} = P^{㊀}_{k} + P_{k}^{㊀}H_k^\top(H_kP_k^{㊀}H_k^\top+R_k)^{-1}(P_K^{㊀}H^\top_k)^\top Pxx=Pk+PkHk(HkPkHk+Rk)1(PKHk) = P k ㊀ + P k ㊀ H k ⊤ ( H k P k ㊀ H k ⊤ + R k ) − 1 ( H k P K ㊀ ) = P^{㊀}_{k} + P_{k}^{㊀}H_k^\top(H_kP_k^{㊀}H_k^\top+R_k)^{-1}(H_kP_K^{㊀}) =Pk+PkHk(HkPkHk+Rk)1(HkPK)

这些本质上是卡尔曼滤波器(或线性MMSE)更新方程。


5. Feature Triangulation

5.1 3D 笛卡尔三角测量

  我们希望创建一个可求解的线性系统,该系统可以让我们初步估计特征的 3D 笛卡尔位置。为此,我们将看到特征的所有位姿都视为已知量。这些特征可以在我们任意选择的相机坐标系 { A } \{A\} { A} 下进行三角化。假设特征点 p f p_f pf 可以在给定的坐标系 { A } \{A\} { A} 下的位姿 1... m 1...m 1...m 下被观测到,我们可以从一些相机位姿 C i , i = 1... m C_i,i=1...m Ci,i=1...m 下得到下列转换。

C i p f = A C i R ( A p f − A p C i ) ^{C_i}p_f = ^{C_i}_AR(^Ap_f-^Ap_{C_i}) Cipf=ACiR(ApfApCi) A p f = A C i R ⊤ C i p f + A p C i ^A{p_f}=^{C_i}_AR^{\top}{^{C_i}p_f}+^Ap_{C_i} Apf=ACiRCipf+ApCi
  在没有噪声的情况下,当前坐标系中的测量值是 方位 C i b ^{C_i}b Cib 及其深度 C i z ^{C_i}z Ciz .因此,我们有以下从当前帧看到的特征的映射:
C i p f = C i z f C i b f = C i z f [ u n v n 1 ] ^{C_i}p_f=^{C_i}z_f^{C_i}b_f = ^{C_i}z_f\left[ \begin{array}{ccc} u_{n} \\ v_{n}\\ 1\\ \end{array} \right] Cipf=CizfCibf=Cizf unvn1

  我们注意到 u n u_n un v n v_n vn表示未畸变的归一化图像坐标。该方位可以通过代入上述公式来转换到固定的坐标系中:
A p f = A C i R ⊤ C i p f + A p C i ^A{p_f}=^{C_i}_AR^{\top}{^{C_i}p_f}+^Ap_{C_i} Apf=ACiRCipf+ApCi = C i z f A b C i → f + A p C i =^{C_i}z_f{^Ab_{C_i\rightarrow f}}+^Ap_{C_i} =CizfAbCif+ApCi

  为了消除估计深度 C i z f ^{C_i}z_f Cizf 额外自由度的需要,我们定义了以下与方位正交的向量 A b C i → f ^Ab_{C_i \rightarrow f} AbCif
A N i = ⌊ A b C i → f × ⌋ = [ 0 − A b C i → f ( 3 ) A b C i → f ( 2 ) A b C i → f ( 3 ) 0 − A b C i → f ( 1 ) − A b C i → f ( 2 ) A b C i → f ( 1 ) 0 ] ^AN_i=\lfloor {^Ab_{C_i \rightarrow f} \times} \rfloor=\left[ \begin{array}{ccc} 0 & -^Ab_{C_i \rightarrow f} (3)& ^Ab_{C_i \rightarrow f}(2)\\ ^Ab_{C_i \rightarrow f} (3) & 0 & -^Ab_{C_i \rightarrow f} (1)\\ -^Ab_{C_i \rightarrow f} (2) & ^Ab_{C_i \rightarrow f} (1) & 0\\ \end{array} \right] ANi=AbCif×= 0AbCif(3)AbCif(2)AbCif(3)0AbCif(1)AbCif(2)AbCif(1)0

  三行都垂直于向量 A b C i → f ^Ab_{C_i \rightarrow f} AbCif,因此 A N i A b C i → f = 0 3 ^AN_i ^Ab_{C_i \rightarrow f}=0_3 ANiAbCif=03。然后,我们可以将变换方程/约束相乘,形成两个仅与未知方程相关的3自由度方程 A p f ^A{p_f} Apf
A N i A p f = A N i C i z f A b C i → f + A N i A p C i ^AN_i^A{p_f}= ^AN_i ^{C_i}{z_f} ^Ab_{C_i \rightarrow f} + ^AN_i ^Ap_{C_i} ANiApf=ANiCizfAbCif+ANiApCi = A N i A p C i =^AN_i^Ap_{C_i} =ANiApCi
通过叠加所有测量值,我们可以:
[ ⋮ A N i ⋮ ] ⏟ A A p f = [ ⋮ A N i A p C i ⋮ ] ⏟ b \underbrace{\left[ \begin{array}{ccc} \vdots \\ ^AN_i\\ \vdots\\ \end{array} \right]}_{A}^A{p_f}=\underbrace{\left[ \begin{array}{ccc} \vdots \\ ^AN_i^Ap_{C_i}\\ \vdots\\ \end{array} \right]}_{b} A ANi pf=b ANiApCi

  由于每个像素测量都提供了两个约束,只要 m > 1 m>1 m>1 ,我们将有足够的约束来对特征进行三角测量。在实践中,特征点的视图越多,三角测量效果越好,因此通常希望从至少五个视图中看到一个特征点。我们可以选择每 A N i ^AN_i ANi 的两行来减少行数,但是通过方形系统,我们可以执行以下“技巧”。
A ⊤ A A p f = A ⊤ b A^\top A ^Ap_f=A^\top b AAApf=Ab ( ∑ i A N i ⊤ A N i ) A p f = ( ∑ i A N i ⊤ A N i A p C i ) (\sum_i{^AN_i^\top {^AN_i}})^Ap_f=(\sum_i{ {^AN^\top _i{^AN_i}{^Ap_{C_i}}}}) (iANiANi)Apf=(iANiANiApCi)

  这是一个 3x3 系统,与原始的 3mx3m 或 2mx2m 系统相比,可以快速解决。我们还检查三角测量特征是否“有效”,并且在相机前方且不太远。上述线性系统和否定系统的 condition number 对错误和极大值数据“敏感”。

5.2 1D 深度三角测量

  我们希望创建一个可求解的线性系统,该系统可以让我们初步估计特征的 1D 深度位置。为此,我们将看到特征的所有姿势与锚架中的方向一起视为已知量。此功能将在我们可以任意选择的某个锚定摄像机帧 { A } \{A\} { A} 中进行三角测量。我们将其定义为锚帧中的归一化图像坐标 [ u n v n 1 ] ⊤ \left[ \begin{array}{ccc} u_{n} &v_{n}&1\\ \end{array} \right]^\top [unvn1]。给定一个锚点姿势 A A A,如果一个特征 p f p_f pf 是可以被位姿 1... m 1...m 1...m 所观测,那么我们可以从任意相机位置 C i , i = 1... m C_i,i=1...m Ci,i=1...m 得到如下转换:
C i p f = A C i R ( A p f − A p C i ) ^{C_i}p_f=^{C_i}_AR(^Ap_f -^Ap_{C_i}) Cipf=ACiR(ApfApCi) A i p f = A C i R ⊤ C i p f + A p C i ^{A_i}p_f=^{C_i}_AR^\top {^{C_i}}p_f +^Ap_{C_i} Aipf=ACiRCipf+ApCi A z f A p f = A C i R ⊤ C i p f + A p C i ^Az_f^{A}p_f=^{C_i}_AR^\top {^{C_i}}p_f +^Ap_{C_i} AzfApf=ACiRCipf+ApCi

  在没有噪声的情况下,当前坐标系中的测量值是方向 C i b ^{C_i}b Cib 及其深度 C i z ^{C_i}z Ciz
C i p f = C i z f C i b f = C i z f [ u n v n 1 ] ^{C_i}p_f=^{C_i}z_f^{C_i}b_f = ^{C_i}z_f\left[ \begin{array}{ccc} u_{n} \\ v_{n}\\ 1\\ \end{array} \right] Cipf=CizfCibf=Cizf unvn1
  我们注意到 u n u_n un v n v_n vn表示未畸变的归一化图像坐标。该方位可以通过代入上述公式来转换到固定的坐标系中:
A z f A b f = A C i R ⊤ C i z f C i b f + A p C i ^A{z_f}^Ab_f=^{C_i}_AR^{\top}{^{C_i}z_f}{^{C_i}b_f} +^Ap_{C_i} AzfAbf=ACiRCizfCibf+ApCi = C i z f A b C i → f + A p C i =^{C_i}z_f{^Ab_{C_i\rightarrow f}}+^Ap_{C_i} =CizfAbCif+ApCi

  为了消除估计深度 C i z f ^{C_i}z_f Cizf 额外自由度的需要,我们定义了以下与方位正交的向量 A b C i → f ^Ab_{C_i \rightarrow f} AbCif
A N i = ⌊ A b C i → f × ⌋ ^AN_i=\lfloor {^Ab_{C_i \rightarrow f} \times} \rfloor ANi=AbCif×
  三行都垂直于向量 A b C i → f ^Ab_{C_i \rightarrow f} AbCif,因此 A N i A b C i → f = 0 3 ^AN_i ^Ab_{C_i \rightarrow f}=0_3 ANiAbCif=03。然后,我们可以将变换方程/约束相乘,形成两个仅与未知方程相关的3自由度方程 A p f ^A{p_f} Apf
( A N i A p f ) A z f = A N i C i z f A b C i → f + A N i A p C i (^AN_i^A{p_f})^Az_f= ^AN_i ^{C_i}{z_f} ^Ab_{C_i \rightarrow f} + ^AN_i ^Ap_{C_i} (ANiApf)Azf=ANiCizfAbCif+ANiApCi = A N i A p C i =^AN_i^Ap_{C_i} =ANiApCi
然后我们可以制定以下系统:
( ∑ i ( A N i A b f ) ⊤ ( A N i A b f ) ) A z f = ( ∑ i ( A N A b f ) i ⊤ A N i A b i ) (\sum_i{(^AN_i^Ab_f)^\top ({^AN_i}}^Ab_f))^Az_f=(\sum_i{ {(^AN^Ab_f)^\top _i{^AN_i}{^Ab_i}}}) (i(ANiAbf)(ANiAbf))Azf=(i(ANAbf)iANiAbi)

  这是一个 1x1 系统,可以用单个标量除法快速求解。我们还检查三角测量特征是否“有效”,并且在相机前方且不太远。完整的特征可以通过以下方式重建:
A p f = A z f A b f ^Ap_f=^Az_f^Ab_f Apf=AzfAbf

5.3 3D Inverse Non-linear Optimization
  获得三角化特征 3D 位置后,将执行非线性最小二乘法来完善此估计值。为了获得良好的数值稳定性,我们对点特征使用逆深度表示,这有助于收敛。我们发现,在大多数情况下,这个问题在室内环境中的 2-3 次迭代内收敛。特征转换可以写成:
C i p f = A C i R ( A p f − A p C i ) ^{C_i}p_f=^{C_i}_AR(^Ap_f -^Ap_{C_i}) Cipf=ACiR(ApfApCi) = A z f A C i R ( [ A x f / A z f A y f / A z f 1 ] − 1 A z f A p C i ) =^Az_f {^{C_i}_AR}\begin{pmatrix} \left[ \begin{array}{ccc} ^Ax_f/^Az_f \\^Ay_f/^Az_f \\ 1\\ \end{array} \right] -{1\over {^Az_f}}^Ap_{C_i}\\ \end{pmatrix} =AzfACiR Axf/AzfAyf/Azf1 Azf1ApCi
   我们定义 u A = A x f / A z f u_A= ^Ax_f/^Az_f uA=Axf/Azf v A = A y f / A z f v_A= ^Ay_f/^Az_f vA=Ayf/Azf ρ A = 1 / A z f \rho_A=1/^Az_f ρA=1/Azf 得到以下测量公式:
h ( u A , v A , ρ A ) = A C i R ( [ u A v A 1 ] − ρ A A p C i ) h(u_A,v_A,\rho_A) = ^{C_i}_AR\begin{pmatrix} \left[ \begin{array}{ccc} u_A\\ v_A \\ 1\\ \end{array} \right] -\rho_A{^Ap_{C_i}}\\ \end{pmatrix} h(uA,vA,ρA)=ACiR uAvA1 ρAApCi

   从相机坐标系看到的特征测量可以重新表述为:
z = [ u i v i ] z=\left[ \begin{array}{ccc} u_i\\ v_i \\ \end{array} \right] z=[uivi] = [ h ( u A , v A , ρ A ) ( 1 ) / h ( u A , v A , ρ A ) ( 3 ) h ( u A , v A , ρ A ) ( 2 ) / h ( u A , v A , ρ A ) ( 3 ) ] =\left[ \begin{array}{ccc} h(u_A,v_A,\rho_A)(1) / h(u_A,v_A,\rho_A)(3)\\ h(u_A,v_A,\rho_A)(2) / h(u_A,v_A,\rho_A)(3)\\ \end{array} \right] =[h(uA,vA,ρA)(1)/h(uA,vA,ρA)(3)h(uA,vA,ρA)(2)/h(uA,vA,ρA)(3)] = h ( u A , v A , ρ A ) =h(u_A,v_A,\rho_A) =h(uA,vA,ρA)

因此,我们可以制定最小二乘法和雅可比:
a r g m i n u A , v A , ρ A ∣ ∣ z − h ( u A , v A , ρ A ) ∣ ∣ 2 \underset{u_A,v_A,\rho_A}{argmin}||z-h(u_A,v_A,\rho_A)||^2 uA,vA,ρAargmin∣∣zh(uA,vA,ρA)2 ∂ h ( u A , v A , ρ A ) h ( u A , v A , ρ A ) = A C i R [ [ 1 0 0 1 0 0 ] − A p C i ] {\partial h(u_A,v_A,\rho_A) \over h(u_A,v_A,\rho_A)} = ^{C_i}_AR\begin{bmatrix} \left[ \begin{array}{ccc} 1 & 0\\ 0&1 \\ 0&0\\ \end{array} \right] -{^Ap_{C_i}}\\ \end{bmatrix} h(uA,vA,ρA)h(uA,vA,ρA)=ACiR 100010 ApCi

   最小二乘问题可以用高斯-牛顿莱文伯格-马夸特算法来解决。


6. Camera Measurement Update

6.1 透视投影(方位)测量模型

   考虑在某个时间 k k k 从相机图像中检测到 3D 特征,其在图像平面上的测量值 u v uv uv(即相应的像素坐标)由下式给出:
z m , k = h ( x k ) + n k = h d ( z n , k , ζ ) + n k = h d ( h p ( C k p f ) , ζ ) + n k = h d ( h p ( h t ( G p f , G G k R , G p C k ) ) , ζ ) + n k = h d ( h p ( h t ( ( h r ( λ , ⋯   ) ) , G G k R , G p C k ) ) , ζ ) + n k z_{m,k} = h(x_k) + n_k \\ =h_d(z_{n,k}, \zeta)+n_k \\ =h_d(h_p(^{C_k}p_f), \zeta)+n_k\\ =h_d(h_p(h_t(^Gp_f,{^{G_k}_GR},^Gp_{C_k})), \zeta)+n_k\\ =h_d(h_p(h_t((h_r(\lambda,\cdots)),{^{G_k}_GR},^Gp_{C_k})), \zeta)+n_k zm,k=h(xk)+nk=hd(zn,k,ζ)+nk=hd(hp(Ckpf),ζ)+nk=hd(hp(ht(Gpf,GGkR,GpCk)),ζ)+nk=hd(hp(ht((hr(λ,)),GGkR,GpCk)),ζ)+nk

   其中 n k n_k nk 是测量噪声,通常假定为零均值白高斯; z n , k z_{n,k} zn,k 是归一化未畸变的 u v uv uv观测值; ζ \zeta ζ是相机的内在参数,如焦距和畸变参数; C k p f ^{C_k}p_f Ckpf是当前相机帧 { C k } \{C_k\} { Ck}中的特征位置; G p f ^{G}p_f Gpf是全局坐标系 { G } \{G\} { G}中的特征位置; { G G k R , G p C k } \{^{G_k}_GR,^Gp_{C_k}\} { GGkR,GpCk}表示全局坐标系中的当前相机姿势(位置和方向)(或相机外参); λ \lambda λ是不同表示(位置除外)的特征参数,例如简单的 xyz 位置或带方位的逆深度。

   在上面的表达式中,我们将测量函数分解为对应于不同操作的多个级联函数,这些函数将状态映射到图像平面上的原始 u v uv uv 测量中。应该注意的是,由于我们将执行具有不同特征表示的内参标定以及外参标定,因此上述相机测量模型是通用的。下一节将给出每个函数的高级说明。

测量函数概述

Function Description
z k = h d ( z n , k , ζ ) z_k=h_d(z_{n,k},\zeta) zk=hd(zn,k,ζ) 将归一化坐标映射到畸变的 UV 坐标的畸变函数
z n , k = h p ( C k p f ) z_{n,k}=h_p(^{C_k}p_f) zn,k=hp(Ckpf) 获取图像中的 3D 点并将其转换为归一化 UV 坐标的投影函数
C k p = h t ( G p f , G C k R , G p C k ) ^{C_k}p_=h_t(^Gp_f,^{C_k}_GR,^Gp_{C_k}) Ckp=ht(Gpf,GCkR,GpCk) 将特征在全局坐标系中的位置转换为当前相机坐标系
G p f = h r ( λ , ⋯   ) ^Gp_f = h_r(\lambda,\cdots) Gpf=hr(λ,) 在全局坐标系中从特征表示转换为 3D 特征

6.2 雅可比计算

   给定上述嵌套函数,我们可以利用链式规则来找到总状态雅可比。由于我们的特征表示函数 h r ( ⋯   ) h_r(\cdots) hr() 也可能取决于状态,即相机姿态,我们需要仔细考虑其附加导数。考虑以下关于雅可比状态 x x x 的测量示例:

∂ z k ∂ x = ∂ h d ( ⋅ ) ∂ z n , k ∂ h p ( ⋅ ) ∂ C k p f ∂ h t ( ⋅ ) ∂ x + ∂ h d ( ⋅ ) ∂ z n , k ∂ h p ( ⋅ ) ∂ C k p f ∂ h t ( ⋅ ) ∂ G P f ∂ h r ( ⋅ ) ∂ x {\partial z_k \over \partial x} = {\partial h_d(\cdot) \over \partial z_{n,k}} {\partial h_p(\cdot) \over \partial ^{C_k}p_f} {\partial h_t(\cdot) \over \partial x} +{\partial h_d(\cdot) \over \partial z_{n,k}} {\partial h_p(\cdot) \over \partial ^{C_k}p_f} {\partial h_t(\cdot) \over \partial ^GP_f}{\partial h_r(\cdot) \over \partial x} xzk=zn,khd()Ckpfhp()xht()+zn,khd()Ckpfhp()GPfht()xhr()

   在全局特征表达(请参阅点特征表达部分)中,第二项将为零,而对于相机表达,则需要计算该项。

6.3 畸变函数

径向畸变
   为了校准相机内参,我们需要知道如何将归一化坐标映射到图像平面上的原始像素坐标。我们首先采用 OpenCV 模型中的径向畸变:

[ u v ] : = z k = h d ( z n , k , ζ ) = [ f x ∗ x + c x f y ∗ y + c y ] \begin{bmatrix} u\\ v \\ \end{bmatrix} :=z_k=h_d(z_{n,k},\zeta) = \begin{bmatrix} f_x*x+c_x \\ f_y*y+c_y \\ \end{bmatrix} [uv]:=zk=hd(zn,k,ζ)=[fxx+cxfyy+cy] x = x n ( 1 + k 1 r 2 + k 2 r 4 ) + 2 p 1 x n y n + p 2 ( r 2 + = 2 x n 2 ) y = y n ( 1 + k 1 r 2 + k 2 r 4 ) + 2 p 2 x n y n + p 1 ( r 2 + = 2 y n 2 ) r 2 = x n 2 + y n 2 x=x_n(1+k_1r^2+k_2r^4)+2p_1x_ny_n+p_2(r^2+=2x^2_n) \\y=y_n(1+k_1r^2+k_2r^4)+2p_2x_ny_n+p_1(r^2+=2y^2_n)\\r^2=x^2_n+y^2_n x=xn(1+k1r2+k2r4)+2p1xnyn+p2(r2+=2xn2)y=yn(1+k1r2+k2r4)+2p2xnyn+p1(r2+=2yn2)r2=xn2+yn2

   其中 z n , k = [ x n y n ] ⊤ z_{n,k} = \begin{bmatrix} x_n & y_n \\ \end{bmatrix}^{\top} zn,k=[xnyn] 是 3D 特征的归一化坐标, u u u v v v 是图像平面上的畸变图像坐标。上述畸变模型涉及以下畸变和相机内(焦距和图像中心)参数,可以在线估算:
ζ = [ f x f y c x c y k 1 k 2 p 1 p 2 ] ⊤ \zeta = \begin{bmatrix} f_x & f_y & c_x & c_y & k_1 & k_2 &p_1 & p_2 \\ \end{bmatrix}^{\top} ζ=[fxfycxcyk1k2p1p2]

   请注意,我们不会像大多数离线校准方法(如 Kalibr)那样估计高阶(即高于四阶)项。要估计这些内参(包括扭曲参数),需要以下这些参数的雅可比矩阵:
∂ h d ( ⋅ ) ∂ ζ = [ x 0 1 0 f x ∗ ( x n r 2 ) f x ∗ ( x n r 4 ) f x ∗ ( 2 x n y n ) f x ∗ ( r 2 + 2 x n 2 ) 0 y 0 1 f y ∗ ( y n r 2 ) f y ∗ ( y n r 4 ) f y ∗ ( r 2 + 2 y n 2 ) f y ∗ ( 2 x n y n ) ] {\partial h_d(\cdot) \over \partial \zeta} = \begin{bmatrix} x & 0 & 1 & 0 & f_x*(x_nr^2) & f_x*(x_nr^4) & f_x*(2x_ny_n) & f_x*(r^2+2x^2_n) \\0 & y & 0 & 1 & f_y*(y_nr^2) & f_y*(y_nr^4) & f_y*(r^2+2y^2_n) & f_y*(2x_ny_n) \end{bmatrix} ζhd()=[x00y1001fx(xnr2)fy(ynr2)fx(xnr4)fy(ynr4)fx(2xnyn)fy(r2+2yn2)fx(r2+2xn2)fy(2xnyn)]

类似地,关于归一化坐标的雅可比坐标可以按如下方式获得:
∂ h d ( ⋅ ) ∂ z n , k {\partial h_d(\cdot) \over \partial z_{n,k}} zn,khd()

鱼眼模型
   由于鱼眼镜头或广角镜头在实践中被广泛使用,我们在这里提供了OpenCV鱼眼镜头中这种畸变模型的数学推导。
[ u v ] : = z k = h d ( z n , k , ζ ) = [ f x ∗ x + c x f y ∗ y + c y ] \begin{bmatrix} u\\ v \\ \end{bmatrix} :=z_k=h_d(z_{n,k},\zeta) = \begin{bmatrix} f_x*x+c_x \\ f_y*y+c_y \\ \end{bmatrix} [uv]:=zk=hd(zn,k,ζ)=[fxx+cxfyy+cy] x = x n r ∗ θ d y = y n r ∗ θ d θ d = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) r 2 = x n 2 + y n 2 θ = a t a n ( r ) x={x_n \over r} * \theta_d \\ y = {y_n \over r} *\theta_d \\ \theta_d = \theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8) \\ r^2 = x^2_n+y^2_n \\ \theta=atan(r) x=rxnθdy=rynθdθd=θ(1+k1θ2+k2θ4+k3θ6+k4θ8)r2=xn2+yn2θ=atan(r)
   其中 z n , k = [ x n y n ] ⊤ z_{n,k} = \begin{bmatrix} x_n & y_n \\ \end{bmatrix}^{\top} zn,k=[xnyn] 是 3D 特征的归一化坐标, u u u v v v 是图像平面上的畸变图像坐标。显然,上述模型中使用了以下畸变内参:
ζ = [ f x f y c x c y k 1 k 2 k 3 k 4 ] ⊤ \zeta = \begin{bmatrix} f_x & f_y & c_x & c_y & k_1 & k_2 &k_3 & k_4 \\ \end{bmatrix}^{\top} ζ=[fxfycxcyk1k2k3k4]

与前面的径向畸变情况类似,校准需要以下这些参数的雅可比矩阵:
∂ h d ( ⋅ ) ∂ ζ = [ x n 0 1 0 f x ∗ ( x n r θ 3 ) f x ∗ ( x n r θ 5 ) ) f x ∗ ( x n r θ 7 ) ) f x ∗ ( x n r θ 9 ) ) 0 y n 0 1 f y ∗ ( y n r θ 3 ) ) f y ∗ ( y n r θ 5 ) ) f y ∗ ( y n r θ 7 ) ) f y ∗ ( y n r θ 9 ) ) ] {\partial h_d(\cdot) \over \partial \zeta} = \begin{bmatrix} x_n & 0 & 1 & 0 & f_x*({x_n \over r}\theta^3) & f_x*({x_n \over r}\theta^5)) & f_x*({x_n \over r}\theta^7)) & f_x*({x_n \over r}\theta^9)) \\0 & y_n & 0 & 1 & f_y*({y_n \over r}\theta^3)) & f_y*({y_n \over r}\theta^5)) & f_y*({y_n \over r}\theta^7)) & f_y*({y_n \over r}\theta^9)) \end{bmatrix} ζhd()=[xn00yn1001fx(rxnθ3)fy(rynθ3))fx(rxnθ5))fy(rynθ5))fx(rxnθ7))fy(rynθ7))fx(rxnθ9))fy(rynθ9))]

类似地,使用微分链式规则,我们可以计算以下关于归一化坐标的雅可比坐标:
∂ h d ( ⋅ ) ∂ z n , k = ∂ u v ∂ x y ∂ x y ∂ x n y n + ∂ u v ∂ x y ∂ x y ∂ r ∂ r ∂ x n y n + ∂ u v ∂ x y ∂ x y ∂ θ d ∂ θ d ∂ θ ∂ θ ∂ r ∂ r ∂ x n y n {\partial h_d(\cdot) \over \partial z_{n,k}} = {\partial uv \over \partial xy} {\partial xy \over \partial x_ny_n} + {\partial uv \over \partial xy} {\partial xy \over \partial r} {\partial r \over \partial x_ny_n} + {\partial uv \over \partial xy} {\partial xy \over \partial \theta_d} {\partial \theta_d \over \partial \theta} {\partial \theta \over \partial r} {\partial r \over \partial x_ny_n} zn,khd()=xyuvxnynxy+xyuvrxyxnynr+xyuvθdxyθθdrθxnynr

6.4 Perspective Projection Function

   标准针孔相机模型用于将相机帧中的 3D 点投影到归一化图像平面(具有单位深度):
z n , k = h p ( C k p f ) = [ C x / C z C y / C z ] z_{n,k} = h_p(^{C_k}p_f) = \begin{bmatrix} ^Cx/{^Cz}\\ ^Cy/{^Cz} \\ \end{bmatrix} zn,k=hp(Ckpf)=[Cx/CzCy/Cz] C k p f = [ C x C y C z ] ^{C_k}p_f= \begin{bmatrix} ^Cx\\ ^Cy \\ {^Cz} \\ \end{bmatrix} Ckpf= CxCyCz
其雅可比矩阵计算如下:
∂ h p ( ⋅ ) ∂ C k p f = [ 1 C z 0 − C x ( C z ) 2 0 1 C z − C y ( C z ) 2 ] {\partial h_p(\cdot) \over \partial ^{C_k}p_f}= \begin{bmatrix} {1 \over {^Cz}} & 0 & {-^Cx\over (^Cz)^2}\\ 0 & {1 \over {^Cz}} & {-^Cy \over (^Cz)^2} \\ \end{bmatrix} Ckpfhp()=[Cz100Cz1(Cz)2Cx(Cz)2Cy]

6.5 欧几里得变换

   我们采用6DOF刚体欧几里得变换,根据当前全局相机姿势将全局坐标系 { G } \{G\} { G} 中的3D特征位置转换为当前相机坐标系 { C k } \{C_k\} { Ck}
C k p f = h t ( G p f , G C k R , G p C k ) = G G k R ( G p f − G p C k ) ^{C_k}p_f = h_t(^Gp_f,{^{C_k}_GR},{^Gp_{C_k} }) = {^{G_k}_GR}({^Gp_f -{^Gp_{C_k}}}) Ckpf=ht(Gpf,GCkR,GpCk)=GGkR(GpfGpCk)

请注意,在视觉惯性导航系统中,我们通常将IMU而不是相机状态保留在状态向量中。因此,我们需要使用 不变的IMU相机外参 { I C R , C p I } \{ {^C_IR}, {^Cp_I}\} { ICR,CpI} 进一步变换上述几何结构,如下所示:
G p C k = G p I k + I G R I p C k = G p I k + I G R I p C {^Gp_{C_k}}={^Gp_{I_k}}+{^G_IR^Ip_{C_k}} = {^Gp_{I_k}} +{^G_IR^Ip_C} GpCk=GpIk+IGRIpCk=GpIk+IGRIpC G C k R = I C k R G I k R = I C R G I k R {^{C_k}_GR } = {^{C_k}_IR} {^{I_k}_GR} = {^{C}_IR} {^{I_k}_GR} GCkR=ICkRGIkR=ICRGIkR

猜你喜欢

转载自blog.csdn.net/qq_39266065/article/details/128529736