前言
空间导航方式分为相对导航Relative Navigation 和迷失空间导航Lost in Space ，相对导航需要依靠其他传感器取得的估计位置信息和姿态信息，依据当前传感器获取的信号完成精调的位姿解算。而迷失空间导航不需要任何初始信息，传感器能够从当前取得的传感信号中独立地解算载体的姿态和位置，具有自主性、无累计误差，是深空探测的必要方式（无论哪种方式，深空探测器都必须要有一种迷失空间导航的传感器，以避免在失去遥测信号后也能够独立地判断当前位置，并与地面站建立联系）。
迷失空间导航
迷失空间导航除了需要当前获取的传感器信号以外，还需要预先存储先验信息，这些先验信息如同人的记忆一样，提供了与导航位置相关的检索目录，传感器通过搜索当前信号与检索目录中信息的匹配程度，获取最可能的目录索引，由此可解算出当前载体最可能的位姿信号，从而完成导航定位。这其中，当前最稳健的导航方式是地形匹配导航，而不变量导航是地形匹配导航中的一种。
陨石坑不变量导航
地外天体最常用、最显著的特征就是其遍布全球的陨石坑，与星表一样，陨石坑可通过预先发送的遥测卫星完成全球陨石坑表绘制，其分辨率可达到米级。陨石坑几何特征明显，能够用简单的数学表示形式存储其特征（通常是圆或者椭圆），减轻了处理器的压力和导航表的存储量。在靠近天体表面时，陨石坑可近似视为平面特征，该特征可通过计算平面的射影不变量完成目录检索，该不变量的定义如下式：
I i j = t r ( C i − 1 C j ) , I j i = t r ( C j − 1 C i ) I_{ij}=\mathrm{tr}(\mathbf{C}_i^{-1}\mathbf{C}_j),I_{ji}=\mathrm{tr}({\mathbf{C}_j^{-1}\mathbf{C}_i})
I ij = tr ( C i − 1 C j ) , I ji = tr ( C j − 1 C i )
式中，C 1 , C 2 \mathbf{C}_1,\mathbf{C}_2 C 1 , C 2 det C 1 = det C 2 = 1 \det{\mathbf{C}_1}=\det{\mathbf{C}_2}=1 det C 1 = det C 2 = 1 [1] 。
陨石坑数学建模
平面上的陨石坑数学描述好建模，那么球面上呢？看上去稍微有点复杂，实际上一点也不简单。
通常绘制的陨石坑目录中，只记录坑中心经纬度坐标、坑的半径，最多最多再记录一个深度信息，不会再有其他的内容了（有也用不上）。
如果要保证记录的目录信息尽量少，同时又保证陨石坑的关键数学模型，那么可以用球面与平面的二维交线作为陨石坑的数学模型，如图所示：
上图说明了两个分别由平面π 1 , π 2 \bm{\pi}_1,\bm{\pi}_2 π 1 , π 2 C 1 , C 2 \mathbf{C}_1,\mathbf{C}_2 C 1 , C 2 l 12 \mathbf{l}_{12} l 12
符号声明
符号名
符号表示
备注
Q \mathbf{Q} Q 陨石坑的三维球方程（矩阵形式）
C \mathbf{C} C 陨石坑的二维平面方程（矩阵形式）
X ~ \mathbf{\tilde{X}} X ~ 非齐次的空间三维点
x ~ \mathbf{\tilde{x}} x ~ 非齐次的平面二维点
X \mathbf{X} X 齐次的空间三维点
x \mathbf{x} x 齐次的平面二维点
l \mathbf{l} l 齐次的二维平面直线
花体符号
对应几何物体点的全集
不共面陨石坑对
不共面的陨石坑对是天体表面常见的陨石坑分布情况。此时不能使用同一个平面计算其不变量，同时，如果将两个陨石坑投影至同一个平面强行制造共面不变量，也将因为投影过程产生的投影畸变使平面不变量数值与真实值不同，因此不能直接套用平面不变量，此时，考虑如上图 所表示的两椭圆所在平面的交线l 12 \mathbf{l}_{12} l 12 C 1 \mathbf{C}_1 C 1 C 2 \mathbf{C}_2 C 2
平面内的二次曲线与两条不平行的直线可形成一个不变量
论文[1] 的着力点就在于此。首先考虑如何利用平面信息求解该空间直线L 12 \mathcal{L}_{12} L 12
C 1 = Q ∩ P 1 , C 2 = Q ∩ P 2 \mathcal{C_1=Q\cap P_1,C_2=Q\cap P_2}
C 1 = Q ∩ P 1 , C 2 = Q ∩ P 2
而直线L 12 = P 1 ∩ P 2 \mathcal{L_{12}=P_1\cap P_2} L 12 = P 1 ∩ P 2
C 1 ∩ C 2 = ( Q ∩ P 1 ) ∩ ( Q ∩ P 2 ) = Q ∩ ( P 1 ∩ P 2 ) = Q ∩ L 12 \mathcal{C_1\cap C_2=(Q\cap P_1)\cap (Q\cap P_2)=Q\cap(P_1\cap P_2)=Q\cap L_{12}}
C 1 ∩ C 2 = ( Q ∩ P 1 ) ∩ ( Q ∩ P 2 ) = Q ∩ ( P 1 ∩ P 2 ) = Q ∩ L 12
因此，空间直线L 12 \mathcal{L}_{12} L 12 C 1 \mathcal{C_1} C 1 C 2 \mathcal{C_2} C 2 l 12 \mathbf{l}_{12} l 12
为求解图像平面上的两个二次曲线交点，联立两个二次曲线的方程：
{ x ⊤ C 1 x = 0 x ⊤ C 2 x = 0 ⟹ x ⊤ ( λ C 1 + C 2 ) x = 0 (1) \left\{
\begin{aligned}
\mathbf{x}^\top\mathbf{C}_1\mathbf{x}&=0\\
\mathbf{x}^\top\mathbf{C}_2\mathbf{x}&=0
\end{aligned}
\right.\implies\mathbf{x}^\top(\lambda\mathbf{C}_1+\mathbf{C}_2)\mathbf{x}=0\tag{1}
{ x ⊤ C 1 x x ⊤ C 2 x = 0 = 0 ⟹ x ⊤ ( λ C 1 + C 2 ) x = 0 ( 1 )
方程(1) 必须是退化的，因为两个二次曲线的交点有且仅有四个（包括复交点）。因此如果将λ C 1 + C 2 \lambda\mathbf{C}_1+\mathbf{C}_2 λ C 1 + C 2
det ( λ C 1 + C 2 ) = 0 ⟹ − C 1 − 1 C 2 − λ I = 0 \begin{aligned}
\det(\lambda\mathbf{C}_1+\mathbf{C}_2)&=0\\
\implies -\mathbf{C}_1^{-1}\mathbf{C}_2-\lambda\mathbf{I}&=0
\end{aligned}
det ( λ C 1 + C 2 ) ⟹ − C 1 − 1 C 2 − λ I = 0 = 0
因此λ \lambda λ − C 1 − 1 C 2 -\mathbf{C}_1^{-1}\mathbf{C}_2 − C 1 − 1 C 2 三个 特征值，其中仅有一个特征值符合要求，下面讨论如何选出该特征值。
求解平面交线
两个不平行的平面有且仅有一条交线，该交线是实的。因此其齐次坐标也应当是实数，对应至方程(1) 的解，应当提供这样的两个实交点的特征值λ 0 \lambda_0 λ 0 B 12 = λ 0 C 1 + C 2 \mathbf{B}_{12}=\lambda_0\mathbf{C}_1+\mathbf{C}_2 B 12 = λ 0 C 1 + C 2 R ( B 12 ) = 2 \mathrm{R}(\mathbf{B}_{12})=2 R ( B 12 ) = 2 B 12 \mathbf{B}_{12} B 12
B 12 = g h ⊤ + h g ⊤ \mathbf{B}_{12} = \mathbf{gh^\top+hg^\top}
B 12 = g h ⊤ + h g ⊤
这里，需要求出两条退化直线g , h \mathbf{g,h} g , h L 12 \mathbf{L}_{12} L 12 l 12 \mathbf{l}_{12} l 12 小结论 ：
三阶二秩矩阵的伴随矩阵可简单地写为：B 12 ∗ = − ( g × h ) ( g × h ) ⊤ = − z z ⊤ \mathbf{B}_{12}^*=-\mathbf{(g\times h)(g\times h)^\top}=-\mathbf{zz^\top} B 12 ∗ = − ( g × h ) ( g × h ) ⊤ = − z z ⊤
因此可知B 12 ∗ \mathbf{B}^*_{12} B 12 ∗ z \mathbf{z} z
B 12 ∗ ≜ [ b 1 ∗ , b 2 ∗ , b 3 ∗ ] ≜ [ b 11 ∗ b 12 ∗ b 11 ∗ b 21 ∗ b 22 ∗ b 23 ∗ b 31 ∗ b 32 ∗ b 33 ∗ ] = − z z ⊤ ⟹ z = − b k ∗ − b k k ∗ \mathbf{B}^*_{12}\triangleq[\mathbf{b}_1^*,\mathbf{b}_2^*,\mathbf{b}_3^*]\triangleq\begin{bmatrix}b_{11}^*&b_{12}^*&b_{11}^*\\b_{21}^*&b_{22}^*&b_{23}^*\\b_{31}^*&b_{32}^*&b_{33}^*\end{bmatrix}=-\mathbf{zz^\top}\implies\mathbf{z}=\frac{-\mathbf{b}_k^*}{\sqrt{-b_{kk}^*}}
B 12 ∗ ≜ [ b 1 ∗ , b 2 ∗ , b 3 ∗ ] ≜ b 11 ∗ b 21 ∗ b 31 ∗ b 12 ∗ b 22 ∗ b 32 ∗ b 11 ∗ b 23 ∗ b 33 ∗ = − z z ⊤ ⟹ z = − b kk ∗ − b k ∗
作者提到，可选取主对角线最大的元素计算以取得更好的数值稳定性[1] z \mathbf{z} z z \mathbf{z} z
B 12 + [ z ] × = 2 g h ⊤ \mathbf{B}_{12}+[\mathbf{z}]_\times=2\mathbf{gh}^\top
B 12 + [ z ] × = 2 gh ⊤
其中任意的非零行是h \mathbf{h} h g \mathbf{g} g l 12 \mathbf{l}_{12} l 12 C 1 , C 2 \mathbf{C}_1,\mathbf{C}_2 C 1 , C 2 x 1 , x 2 \mathbf{x}_1,\mathbf{x}_2 x 1 , x 2 l 12 \mathbf{l}_{12} l 12
l 12 ⊤ x 1 l 12 ⊤ x 2 < 0 \mathbf{l}^\top_{12}\mathbf{x}_1\mathbf{l}^\top_{12}\mathbf{x}_2<0
l 12 ⊤ x 1 l 12 ⊤ x 2 < 0
满足该要求的解为所求直线射影l 12 \mathbf{l}_{12} l 12
因此，该特征值λ 0 \lambda_0 λ 0 λ 0 \lambda_0 λ 0 l 12 \mathbf{l}_{12} l 12
求解不变量
平面上的二次曲线C \mathbf{C} C x , y \mathbf{x,y} x , y
I = ( x ⊤ C y ) 2 ( x ⊤ C x ) ( y ⊤ C y ) I=\frac{(\mathbf{x^\top C y})^2}{(\mathbf{x^\top C x})(\mathbf{y^\top C y})}
I = ( x ⊤ Cx ) ( y ⊤ Cy ) ( x ⊤ Cy ) 2
同理，由于二维平面上的直线与点存在对偶关系，下式亦为射影不变量：
J = ( m ⊤ C ∗ n ) 2 ( m ⊤ C ∗ m ) ( n ⊤ C ∗ n ) J=\frac{(\mathbf{m^\top C^* n})^2}{(\mathbf{m^\top C^* m})(\mathbf{n^\top C^* n})}
J = ( m ⊤ C ∗ m ) ( n ⊤ C ∗ n ) ( m ⊤ C ∗ n ) 2
则根据以上结果，可利用三个不共面的陨石坑C 1 , C 2 , C 3 \mathbf{C}_1,\mathbf{C}_2,\mathbf{C}_3 C 1 , C 2 , C 3 l 12 , l 13 , l 23 \mathbf{l}_{12},\mathbf{l}_{13},\mathbf{l}_{23} l 12 , l 13 , l 23
{ α 1 2 = ( l 12 ⊤ C 1 ∗ l 13 ) 2 ( l 12 ⊤ C 1 ∗ l 12 ) ( l 13 ⊤ C 1 ∗ l 13 ) α 2 2 = ( l 23 ⊤ C 2 ∗ l 12 ) 2 ( l 23 ⊤ C 2 ∗ l 23 ) ( l 12 ⊤ C 2 ∗ l 12 ) α 3 2 = ( l 13 ⊤ C 3 ∗ l 23 ) 2 ( l 13 ⊤ C 3 ∗ l 13 ) ( l 23 ⊤ C 3 ∗ l 23 ) \left\{
\begin{aligned}
\alpha_1^2=\frac{(\mathbf{\mathbf{l}_{12}^\top C^*_1 \mathbf{l}_{13}})^2}{(\mathbf{\mathbf{l}_{12}^\top C^*_1 \mathbf{l}_{12}})(\mathbf{\mathbf{l}_{13}^\top C^*_1 \mathbf{l}_{13}})}\\
\alpha_2^2=\frac{(\mathbf{\mathbf{l}_{23}^\top C^*_2 \mathbf{l}_{12}})^2}{(\mathbf{\mathbf{l}_{23}^\top C^*_2 \mathbf{l}_{23}})(\mathbf{\mathbf{l}_{12}^\top C^*_2 \mathbf{l}_{12}})}\\
\alpha_3^2=\frac{(\mathbf{\mathbf{l}_{13}^\top C^*_3 \mathbf{l}_{23}})^2}{(\mathbf{\mathbf{l}_{13}^\top C^*_3 \mathbf{l}_{13}})(\mathbf{\mathbf{l}_{23}^\top C^*_3 \mathbf{l}_{23}})}
\end{aligned}
\right.
⎩ ⎨ ⎧ α 1 2 = ( l 12 ⊤ C 1 ∗ l 12 ) ( l 13 ⊤ C 1 ∗ l 13 ) ( l 12 ⊤ C 1 ∗ l 13 ) 2 α 2 2 = ( l 23 ⊤ C 2 ∗ l 23 ) ( l 12 ⊤ C 2 ∗ l 12 ) ( l 23 ⊤ C 2 ∗ l 12 ) 2 α 3 2 = ( l 13 ⊤ C 3 ∗ l 13 ) ( l 23 ⊤ C 3 ∗ l 23 ) ( l 13 ⊤ C 3 ∗ l 23 ) 2
作者进一步指出，如果能够保证导航坑对都不直接相交，则以上三组不变量可以均可以用正的、数值幅度增大的复合不变量表示：
{ J 1 = a r c c o s h α 1 = ln ( α 1 + α 1 2 − 1 ) J 2 = a r c c o s h α 2 = ln ( α 2 + α 2 2 − 1 ) J 3 = a r c c o s h α 3 = ln ( α 3 + α 3 2 − 1 ) \left\{
\begin{aligned}
J_1=\mathrm{arccosh}{\alpha_1}=\ln(\alpha_1+\sqrt{\alpha_1^2-1})\\
J_2=\mathrm{arccosh}{\alpha_2}=\ln(\alpha_2+\sqrt{\alpha_2^2-1})\\
J_3=\mathrm{arccosh}{\alpha_3}=\ln(\alpha_3+\sqrt{\alpha_3^2-1})
\end{aligned}
\right.
⎩ ⎨ ⎧ J 1 = arccosh α 1 = ln ( α 1 + α 1 2 − 1 ) J 2 = arccosh α 2 = ln ( α 2 + α 2 2 − 1 ) J 3 = arccosh α 3 = ln ( α 3 + α 3 2 − 1 )
通过选取陨石坑三元组构成导航坑表，检索当前图像中陨石坑三元组在表中三元组的位置，即可完成三元组的定位。
一定程度上，三元组的检索甚至快于二元组，因为它直接就能构成为金字塔，省去了从二元组计算三元组的步骤。
参考文献
