小天体探测有助于揭示太阳系的起源与演化，并为行星防御积累技术和经验，具有重大的科学意义和应用价值^[1-4]。迄今为止，世界各国已经开展了接近20次小天体探测任务^[5-9]，包括飞掠、撞击、环绕探测和采样返回等形式，获得了丰富的科学数据并实现了巨大的技术突破。然而，小天体与地球距离遥远，且其引力小且不规则，这严重制约了探测器安全并高效地执行任务，同时对导航系统的自主性、稳定性和准确性提出了非常严苛的要求。

利用小天体的光学图像与已知的表面陆标进行匹配，并建立2D到3D信息的对应关系^[10]，可以为探测器提供关键的定位信息，从而有望实现精确的自主导航。由此可见，高效、鲁棒的特征/陆标提取与匹配是成功实施未来小天体探测任务的使能技术之一^[11]。由于小天体表面分布着大量的岩石碎片与沟壑等天然地形，可从中提取出丰富的纹理特征，如角点和曲线特征等。角点特征的提取算法包括SIFT（Scale Invariant Feature Transform）^[12]、SURF（Speeded Up Robust Features）^[13]和ORB（Oriented fast and Rotated Brief）^[14]等，通过计算对应特征描述子之间的距离来匹配。此外，崔平远等^[15]利用激光雷达获取的高程图构建了三维地形的几何不变特征，并提出了基于投票的特征匹配策略；邵巍等^[16]提取了图像信息更丰富的曲线特征，采用最近邻距离比率的方法匹配；王光泽等^[17]针对曲线特征提出了基于曲率的匹配算法。另一方面，陨石坑作为一种显著的天然陆标，学者对其检测与匹配方法进行了深入的研究。Bilodeau等^[18]提出了一种基于层次分割的陨石坑检测与匹配方法；冯军华等^[19]利用Canny算子和边缘配对实现陨石坑的检测，并提出了基于投票策略的陨石坑匹配方法；邵巍等^[20]利用深度学习对陨石坑进行智能检测，较好地解决了小天体暗弱环境下陨石坑的漏检问题；Zhu等^[21]利用最小二乘法对陨石坑的边缘进行拟合与定位，并给出了其误差的不确定性；Shao等^[22]利用圆弧支撑域构造了陨石坑的描述子，并提出了一种采用最近邻距离比和欧氏距离约束的陨石坑匹配算法。

然而，小天体附近极端的观测条件对陆标匹配的鲁棒性提出了极高的要求。由于小天体周围不存在大气覆盖，没有光线照射的区域在图像中将呈现出阴影；在不同距离和不同视角条件下，小天体在图像中将呈现出不同的视觉外观；小天体表面可能不存在足够多的陨石坑用于导航，需要研究广义陆标（一般性的表面地形）的鲁棒匹配方法。

“隼鸟号”（ Hayabusa）^[23-24]通过GCP-NAV软件匹配L-map表面陆标来导航，“罗塞塔号”（Rosetta）^[25]和“欧西里斯–雷克斯”（Origins Spectral Interpretation Resource Identification Security Regolith Explorer，OSIRIS-REx）^[26]通过匹配maplet陆标导航。不管是L-map还是maplet，其本质都是小天体上一块局部的高精度地形/反照率图，由立体光度测量（Stereo-Photo Clinometry，SPC）法^[27]利用在不同光照、不同视角条件下拍摄的多幅图像通过迭代计算而得到。在给定的观测和光照条件下，可以渲染出该陆标的视觉外观，并利用归一化互相关（Normalized Cross-Correlation，NCC）算法将其与实际拍摄的图像匹配。该方法可以适应极端的尺度、视角和光照变化，匹配精度高且理论上可应用于任何地形。但计算复杂，难以在星载计算机中实时处理，且匹配结果对相机的位姿估计误差十分敏感。

为在小天体附近的极端环境中为探测器提供准确可靠的视觉信息，本文提出了针对L-map/maplet陆标（以下简称“陆标”）的鲁棒匹配方法。首先，介绍了陆标的光度模型及其在图像中的像素投影模型，并基于像素投影模型对陆标的匹配误差进行了理论分析。然后，基于误差分析对陆标点进行优化，并提出了加权归一化互相关（Weighted Normalized Cross-Correlation，WNCC）算法以对陆标进行鲁棒匹配。最后，利用合成图像序列对比分析了WNCC与NCC算法在系统误差、尺度、视角和太阳相位角变化等条件下的匹配性能。

本文讨论的陆标是由SPC技术生成的局部、高精度的地形/反照率图^[27]，也称为L-map或maplet，由中心（特征）点及其周围的陆标点构成。以陆标的中心点为原点建立东–北–天陆标坐标系$ {\mathcal{F}^{\text{L}}} $，如图1所示。利用从不同角度和光照条件下拍摄的小天体图像，SPC技术可以估计出各陆标点在$ {\mathcal{F}^{\text{L}}} $中的高程及其所在表面的法向量和反照率。

图  1  由SPC技术生成的陆标及其在小天体中的位置

Figure 1.  Landmark generated by SPC technology and its position on asteroid surface

陆标匹配过程分为3步：第1步，在给定的观测（相机与陆标的相对位姿）和光照条件下，通过光度模型渲染出陆标在相机中的视觉外观，本文将其定义为渲染图像$ {{\boldsymbol{I}}_{{\text{ren}}}} $。第2步，根据观测条件通过投影模型估计陆标上各点在实际拍摄图像${\boldsymbol{ I}} $当中的像素坐标，得到映射坐标构成的张量矩阵${\boldsymbol{ M}} $。第3步，在$ {\boldsymbol{I}} $中的一定范围内整体平移$ {\boldsymbol{M}} $得到新的张量矩阵，并利用插值计算每一个陆标点在${\boldsymbol{ I}} $中对应的亮度值，将得到的图像定义为校正图像$ {{\boldsymbol{I}}_{{\text{rec}}}} $；重复执行第3步，直到$ {{\boldsymbol{I}}_{{\text{ren}}}} $与$ {{\boldsymbol{I}}_{{\text{rec}}}} $的NCC值最大，则陆标中心点所对应的像素坐标即为匹配结果。

在真空环境和太阳平行光照射条件下，陆标上一点$ {\boldsymbol{p}} $通过相机被转换为像素亮度值的过程可以通过光度模型来描述。定义$ {\boldsymbol{p}} $点所在的表面法向量为${\boldsymbol{ n}} $，太阳光入射方向为$ {{\boldsymbol{i}}_{\text{s}}} $、反射方向$ {e_{\text{s}}} $（即观测相机的方向），如图2所示。用$ {i_{\text{s}}} $表示$ {{\boldsymbol{i}}_{\text{s}}} $与${\boldsymbol{ n}} $的夹角（即入射角），$ {e_{\text{s}}} $表示$ {{\boldsymbol{e}}_{\text{s}}} $与$ {\boldsymbol{n}} $的夹角（即反射角），α为${{\boldsymbol{e}}_{\text{s}}}$与$ {{\boldsymbol{i}}_{\text{s}}} $的夹角（即太阳相位角），则点$ {\boldsymbol{p}} $在图像中的像素亮度$ {{\boldsymbol{I}}_{\text{p}}} $可建模为

图  2  光度模型中的向量及角度定义

Figure 2.  Definitions of vectors and angles in photometry

其中：$ \Lambda \left( {{I_{\text{s}}},{t_{\text{i}}},{K_{\text{C}}}} \right) $是太阳光强度$ {I_{\text{s}}} $、相机曝光时间$ {t_{\text{i}}} $和相机参数$ {K_{\text{C}}} $（$ {K_{\text{C}}} $表示将入射光的强度转换为相机CCD 像素亮度大小的比例因子）的函数；$ \mathcal{A} $为${\boldsymbol{p}}$点的表面反照率；$ B $为图像背景亮度；$ R\left( {\alpha ,{i_{\text{s}}},{e_{\text{s}}}} \right) $是与图2中所示角度相关的光线反射函数，具体可参考文献[27]。根据式（1）可以计算出所有陆标点的像素亮度值，得到的图像即称为渲染图像${{\boldsymbol{I}}_{{\text{ren}}}}$。

在小天体质心固连坐标系$ {\mathcal{F}^{\text{G}}} $下，假设相机的空间位置$ r \in {\mathbb{R}^3} $，其姿态可以利用李代数表示为${\boldsymbol{\xi}} =$$ {\left[ {\begin{array}{*{20}{c}} {{\xi _x}}&{{\xi _y}}&{{\xi _z}} \end{array}} \right]^{\text{T}}} \in \mathfrak{s}\mathfrak{o}\left( 3 \right) $。从相机坐标系$ {\mathcal{F}^{\text{C}}} $到$ {\mathcal{F}^{\text{G}}} $的旋转矩阵可表示为${{\boldsymbol{R}}_{{\text{c2g}}}} = \exp \left({{\boldsymbol{\xi}}^{\wedge}} \right)$，其中${{\boldsymbol{\xi}}^{\wedge}}$为$ \xi $中各元素所构成的反对称矩阵。

陆标的像素投影模型如图3所示。假设陆标的中心点在$ {\mathcal{F}^{\text{G}}} $中的位置为$ L \in {\mathbb{R}^3} $，从$ {\mathcal{F}^{\text{L}}} $到$ {\mathcal{F}^{\text{G}}} $的旋转矩阵为${{\boldsymbol{R}}_{{\text{l2g}}}}$，则陆标中心点在$ {\mathcal{F}^{\text{C}}} $中的位置${\boldsymbol{V}}$为

图  3  陆标的像素投影模型

Figure 3.  Pixel projection model of landmark

令${{\boldsymbol{L}}_{i}} \in {\mathbb{R}^3}$表示任意陆标点在$ {\mathcal{F}^{\text{L}}} $中的位置，其在$ {\mathcal{F}^{\text{C}}} $中的坐标${{\boldsymbol{V}}_i}$为

根据式（3）和式（4）可计算出陆标中心点和各陆标点在像平面内的像素坐标${\boldsymbol{u}} \in {\mathbb{R}^2}$，${{\boldsymbol{u}}_i} \in {\mathbb{R}^2}$

其中：$\bar {\boldsymbol{u}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}^{\rm T}}}&1 \end{array}} \right]^{\rm T}}$，${\bar {\boldsymbol{u}}_{\text{i}}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{u}}_{\text{i}}^{\text{T}}}&1 \end{array}} \right]^{\text{T}}}$分别为${\boldsymbol{u}}$和${{\boldsymbol{u}}_{i}}$的齐次坐标；$ {s_1} $和$ {s_2} $为尺度因子；${\boldsymbol{K}}$为相机的内参矩阵。

其中：$ f $为相机的焦距；$ {u_x} $和$ {u_y} $为相机主点在图像坐标系中的像素偏移量。

则${\boldsymbol{L}}$和${{\boldsymbol{L}}_i}$在${\boldsymbol{I}}$中所对应的像素坐标构成了张量矩阵${\boldsymbol{M}}$，通过像素插值可以得到校正图像${{\boldsymbol{I}}_{{\text{rec}}}}$。

根据陆标匹配过程可知，匹配误差主要来源于光度和像素投影模型误差。然而，光度模型是实际光线反射过程的近似，无法对复杂的多次反射和散射进行准确的建模，故难以对其进行误差分析。本节主要基于陆标的像素投影模型分析陆标匹配误差产生的机制，并为后续陆标点的优化及陆标鲁棒匹配算法的设计提供依据。

考虑到陆标的位置及相机的位姿先验信息均存在噪声，利用上述陆标投影模型将会导致陆标在图像中的投影区域（即${\boldsymbol{M}}$）发生变形，从而造成匹配误差，匹配误差产生的机制如图4所示。由于陆标匹配的目的是在图像${\boldsymbol{I}}$中寻找陆标中心点${\boldsymbol{L}}$所对应的像素坐标，故可以将${\boldsymbol{L}}$在图像中的投影作为参考点，利用陆标点${{\boldsymbol{L}}_i}$与${\boldsymbol{L}}$之间的像素坐标偏移误差来描述陆标在图像中投影区域的变形程度。

图  4  陆标匹配误差产生的机制

Figure 4.  Mechanism of landmark matching error

此后，用$\hat {\boldsymbol{*}}$指代向量*的估计值。假设${{\boldsymbol{L}}_i}$与${\boldsymbol{L}}$在图像中像素坐标的偏移${{\boldsymbol{d}}_i} = {{\boldsymbol{u}}_i} - {\boldsymbol{u }}$，为描述${\boldsymbol{L}}$在$ {\mathcal{F}^{\text{G}}} $中的位置噪声$\Delta {\boldsymbol{L}}\left( {\hat {\boldsymbol{L}} = {\boldsymbol{L}} + \Delta{\boldsymbol{ L}}} \right)$、${{\boldsymbol{L}}_i}$在$ {\mathcal{F}^{\text{L}}} $中的相对位置噪声$\Delta {{\boldsymbol{L}}_i}\left( {{{\hat {\boldsymbol{L}}}_i} = {{\boldsymbol{L}}_i} + \Delta {{\boldsymbol{L}}_i}} \right)$及相机在$ {\mathcal{F}^{\text{G}}} $中的位置估计噪声$\Delta {\boldsymbol{r}}\left( {\hat {\boldsymbol{r}} = {\boldsymbol{r}} + \Delta {\boldsymbol{r}}} \right)$和姿态估计噪声$\Delta {\boldsymbol{\xi}} \left( {\hat {\boldsymbol{\xi}} = {\boldsymbol{\xi}} + \Delta {\boldsymbol{\xi}} } \right)$对${{\boldsymbol{d}}_i}$的影响，推导了${{\boldsymbol{d}}_i}$相对于$\Delta {\boldsymbol{L}}$、$\Delta {{\boldsymbol{L}}_i}$、$\Delta {\boldsymbol{r}}$和$\Delta {\boldsymbol{\xi}}$的偏导数，得到雅可比矩阵${{\boldsymbol{J}}_i}$为

其中：

令上述误差项所引起的${{\boldsymbol{d}}_i}$的增量为$\Delta {{\boldsymbol{d}}_i}$，即${{\boldsymbol{L}}_i}$相对于陆标中心点${\boldsymbol{L}}$在图像中像素投影的变形量，则有

式（9）描述了陆标的位置误差和相机的位姿误差如何影响陆标在图像中的变形程度。例如，在完全相同的误差条件下，当相机的位置不同时，${{\boldsymbol{L}}_i}$相对于${\boldsymbol{L}}$的变形程度如图5所示。其中，图5（a）的相机在$ {\mathcal{F}^{\text{G}}} $中的位置为$[ {\begin{array}{*{20}{c}} { - 56.7}&{ - 89.4}&{151.1} \end{array}} ]{\;{\text{m}}} $，图5（b）相机的位置为$ [ {\begin{array}{*{20}{c}} { - 63.3}&{ - 76.8}&{135.7} \end{array}} ]\;{\text{m}} $。可以看出，距中心点越远，陆标点的变形程度受相机和陆标噪声的影响越大；此外，不同的相机状态也会影响其误差传播。这些发生严重变形的点对应图像${\boldsymbol{I}}$中错误的像素坐标，使得到的校正图像${{\boldsymbol{I}}_{{\text{rec}}}}$与渲染图象${{\boldsymbol{I}}_{{\text{ren}}}}$在外观上存在很大的不同，故导致了匹配误差。利用带有误差的参数将陆标投影到图像平面上，得到陆标的实际像素区域与真实像素区域不重合。

图  5  陆标上各点相对于中心点的像素变形程度

Figure 5.  Pixel deformation of landmark

根据以上分析可知，当相机的状态或陆标的位置包含较大噪声时，与陆标中心点相距较远的陆标点在图像中的位置估计结果达到了几个甚至几十个像素误差，利用这些点不但不能提升匹配精度，反而会影响匹配结果；此外，陆标的像素变形情况与相机的状态存在紧密联系，即在相同的噪声情况下，不同相机状态将导致不同的误差分布。因此，动态地选取陆标表面上一些变形程度较小的陆标点生成${{\boldsymbol{I}}_{{\text{ren}}}}$和${{\boldsymbol{I}}_{{\text{rec}}}}$并进行匹配将会提高匹配精度，同时根据每个点的变形程度动态地分配不同的权值也能降低陆标变形对匹配结果的影响，从而提升鲁棒性。

假设$\Delta {\boldsymbol{L}}$、$\Delta {{\boldsymbol{L}}_{i}}$、$\Delta {\boldsymbol{r}}$和$\Delta {\boldsymbol{\xi}}$的各个分量均为独立同分布的高斯白噪声，其标准差已知，分别为${\sigma _{L}}$、${\sigma _{{L,i}}}$、${\sigma _{r}}$和$ {\sigma _\xi } $。为了量化描述陆标点的变形程度$\Delta {{\boldsymbol{d}}_i}$在当前相机状态下受$\Delta {\boldsymbol{L}}$、$\Delta {{\boldsymbol{L}}_i}$、$\Delta {\boldsymbol{r}}$和$ \Delta \xi $的影响，推导了其协方差并定义了描述其变形的度量因子。根据式（9）可以计算出$\Delta {{\boldsymbol{d}}_i}$的协方差

其中：${{\boldsymbol{I}}_3}$表示3阶单位矩阵，${{\boldsymbol{0}}_3}$表示大小为$ 3 \times 3 $的零矩阵。定义度量因子$ {\delta _i} $（单位：pixel）为${\bf{E}}\left\{ {\Delta {{\boldsymbol{d}}_i}} \right\}$的迹的平方根，记作

则$ {\delta _i} $表征了在给定的噪声条件下，通过像素投影模型，陆标点${{\boldsymbol{L}}_i}$相对于中心点L在图像中像素的变形程度，其值越大，说明$ {{\boldsymbol{L}}_i} $越容易受到上述误差项的影响，所以利用$ {{\boldsymbol{L}}_i} $来进行陆标匹配将会造成更大的匹配误差。

为降低计算复杂度并提高陆标的匹配精度和鲁棒性，设计了陆标点优化选取方法以及WNCC陆标鲁棒匹配方法。

“罗塞塔号”（Rosetta）任务使用的maplet陆标大小为$99 \times 99$，共由9 801个点组成^[25]。理论上利用所有陆标点的信息可以提高匹配结果的精度和稳定性，然而在特定的观测条件及噪声条件下，其中一些点可能存在较大的像素投影误差，不适合用于陆标匹配。为了降低计算复杂度并消除这些点对匹配结果的影响，根据度量因子$ {\delta _i} $对陆标点进行优化选取得到集合${{\boldsymbol{S}}_L}$为

其中：$ {T_\delta } $（单位：pixel）为度量因子的阈值。

然而，由于地形遮挡，太阳入射光线无法照射的区域将在图像当中呈现出阴影，而其反射光线进入相机视场的陆标点将不可见，如图6所示，其中红色表示在相机中可见的光照区域，黑色表示阴影，灰色表示不可见区域。黑色和灰色区域所覆盖的陆标点都应该从${{\boldsymbol{S}}_L}$中剔除。

图  6  地形遮挡造成的阴影与不可见区域

Figure 6.  Shadows and invisible areas caused by terrain occlusion

为判断陆标点在相机当中是否可见，可以计算该点与相机的连线之间是否与陆标网格模型存在其他交点，例如图6中线段$\overline {{{\boldsymbol{L}}_i}{\boldsymbol{r}}}$与陆标之间存在交点${{\boldsymbol{L}}_j}$，故$ {{\boldsymbol{L}}_i} $不可见，而$ {{\boldsymbol{L}}_j} $可见。同理，若陆标点与太阳之间的连线与陆标网格模型存在其他交点，则该点为阴影，如${{\boldsymbol{L}}_{\rm p}}$。

此外，若${{\boldsymbol{S}}_L}$中任意两个陆标点在图像中的像素投影距离太近，那么这两点所对应的像素亮度值相近，不具备良好的区分度。可以设置距离阈值${T_{\rm dis}}$，当两点之间的像素距离小于$ {T_{{\text{dis}}}} $时，则剔除其中度量因子较大的一个。

为提高陆标匹配的鲁棒性，根据每一个被选取的点${{\boldsymbol{L}}_i} \in {{\boldsymbol{S}}_L}$所对应的度量因子$ {\delta _i} $，定义其在陆标匹配中的权值$ {w_i} $（置信度）

其中：s为尺度因子。

假设选取的陆标点数量为$n = {\rm size}\left( {{{\boldsymbol{S}}_L}} \right)$，则由式（1）可以渲染出$ {{\boldsymbol{S}}_L} $中每一点$ {L_i} $的灰度值$ I_{{\text{ren}}}^i $，得到向量形式的渲染图像$ {{\boldsymbol I}_{{\text{ren}}}} = \left[ {\begin{array}{*{20}{c}} {I_{{\text{ren}}}^{\text{1}}}&{I_{{\text{ren}}}^{\text{2}}}& \cdots &{I_{{\text{ren}}}^n} \end{array}} \right] $，由式（4）可以计算出${{\boldsymbol{L}}_i}$在图像中的像素坐标${{\boldsymbol{u}}_i}$，得到映射坐标张量${\boldsymbol M} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_1}}&{{{\boldsymbol{u}}_2}}& \cdots &{{{\boldsymbol{u}}_n}} \end{array}} \right]$。利用${\boldsymbol M} $通过像素插值可以计算${{\boldsymbol{u}}_i}$在图像${\boldsymbol{I}}$中对应的亮度值$I_{\rm rec}^i$，得到向量形式的校正图像$ {{\boldsymbol{I}}_{\rm rec}} = \left[ {\begin{array}{*{20}{c}} {I_{\rm rec}^1}&{I_{\rm rec}^2}& \cdots &{I_{\rm rec}^n} \end{array}} \right] $，根据式（13）可以计算出权值向量${\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}} {{w_1}}&{{w_2}}& \cdots &{{w_n}} \end{array}} \right]$。则${{\boldsymbol I}_{\rm ren}}$和${\boldsymbol{I}_{\rm rec}}$之间的加权归一化互相关值$ {C_w} $（范围：[−1, 1]）可以通过下式计算得到

其中：

运算符∑表示对向量中的所有元素求和，⊙表示向量中对应元素相乘。在一定范围内平移映射坐标张量$ {\boldsymbol{M}} $并利用像素插值得到新的$ {{\boldsymbol{I}}_{{\text{rec}}}} $，再计算其与$ {{\boldsymbol{I}}_{{\text{ren}}}} $之间的归一化互相关值。不断重复该过程，记录其中最大$ {C_w} $值所对应的像素位置作为陆标在图像$ {\boldsymbol{I}} $中的匹配结果。

为防止在全图中搜索、匹配陆标, 以提高算法效率，需要对搜索范围进行动态估计。根据式（3）推导陆标中心$ {\boldsymbol{L}} $在图像$ {\boldsymbol{I}} $中的像素坐标$ {\boldsymbol{u}} $相对于$ \Delta {\boldsymbol{L}} $、$ \Delta{\boldsymbol{ r}} $和$\Delta {\boldsymbol{\xi }}$的偏导数，得到雅可比矩阵$ {{\boldsymbol{J}}_L} $

则$ {\boldsymbol{u}} $在图像$ {\boldsymbol{I}} $中像素估计误差$ \Delta {\boldsymbol{u}} $（$ \hat {\boldsymbol{u}} = {\boldsymbol{u}} + \Delta{\boldsymbol{ u}} $）的协方差可估计为

定义搜索区域的半径为${R_{\rm s}}$（pixels）

其中：$ {R_0} $为最小搜索半径，其目的是为了防止${R_{\rm s}}$过小导致搜索区域无法覆盖陆标在$ {\boldsymbol{I}} $中的真实位置。

为了在图像中得到亚像素级的陆标定位，通过曲面拟合计算相关值$ {C_w} $的极值所对应的像素坐标，此处不进行详细讨论。WNCC陆标鲁棒匹配算法的流程如图7所示。

图  7  WNCC陆标匹配算法流程

Figure 7.  Flowchart of WNCC landmark matching algorithm

为验证WNCC 陆标匹配算法的性能，开发了小天体视觉仿真系统（Visual Simulation System, VSS）合成了如图1所示的陆标（其在$ {\mathcal{F}^{\text{G}}} $中的位置为$ {\boldsymbol L} = {[ {\begin{array}{*{20}{c}} { - 38.3}&{ - 60.3}&{101.1} \end{array}} ]^{\text{T}}}{\;\text{m}} $，网格大小为$ 99 \times 99 $，分辨率：$ 0.3\;{\text{m}} $）在不同距离、不同视角和不同光照等极端条件下的序列图像，并在完全相同的噪声条件下比较了WNCC算法和目前小天体探测任务^[7-10]使用的NCC算法的匹配性能。

设计了7组仿真算例，分别讨论了陆标点位置误差、相机与陆标之间的相对位置误差、姿态误差、陆标点数量变化、图像尺度变化、相机视角变化和太阳相位角变化对匹配结果的影响。对于前4组仿真算例，设置光照方向为$[ {\begin{array}{*{20}{c}} 0&{ - 1}&0 \end{array}}] $，相机在$ {\mathcal{F}^{\text{L}}} $下的位置为$ [ {\begin{array}{*{20}{c}} 0 &0 &{200} \end{array}} ]\;{\;\text{m}} $，姿态为$[ {\begin{array}{*{20}{c}} {180^\circ } &{0^\circ } &{0^\circ } \end{array}} ] $（3-1-2旋转序列的欧拉角），且在空间当中保持不变，相机拍摄的图像如图8所示。其中，前3组算例设置WNCC选取$ n = 500 $个陆标点进行陆标匹配。对于后3组仿真算例，设置标称系统噪声参数如表1所示。每组算例均执行1 000 次蒙特卡洛仿真，得到的统计结果以陆标像素均方根误差（Root Mean Squared Error，RMSE）的形式呈现。值得注意的是，尽管式（1）所描述的光度模型与实际情况存在一定误差，但其对陆标匹配结果的影响非常小，可忽略不计。

图  8  相机拍摄的图像（与陆标相距200 m）

Figure 8.  Image captured by the camera (200 m from the landmark)

首先，考虑陆标点$ {{\boldsymbol{L}}_i} $在$ {\mathcal{F}^{\text{L}}} $中的位置误差对匹配结果的影响。令$ {\sigma _{L,i}} $从0.05 m开始且以0.05 m的间隔均匀增加到1.05 m，$ {\sigma _L} = 0\;{\text{m}} $、$ {\sigma _r} = 0\;{\text{m}} $且$ {\sigma _\xi } = 0^\circ $。在不同的陆标点位置误差影响下，WNCC和NCC算法的陆标匹配结果如图9所示。

图  9  在不同$ {\sigma _{L,i}} $影响下，WNCC与NCC的匹配结果

Figure 9.  Landmark matching results of WNCC and NCC under different $ {\sigma _{L,i}} $

可以看出，随着$ {\sigma _{L,i}} $的增大，WNCC和NCC算法的匹配精度逐渐降低，但WNCC算法具有更高的精度，对陆标点的位置误差具有更强的稳定性。

考虑到陆标中心点L的位置误差与相机位置误差是相对的，故只需要讨论其中一种情况，此处只分析相机位置误差对匹配结果的影响。令$ {\sigma _r} $由0.1 m以0.2 m的间隔增加到4.1 m，$ {\sigma _L} = 0\;{\text{m}} $、$ {\sigma _{L,i}} = 0\;{\text{m}} $且$ {\sigma _\xi } = 0^\circ $。在不同的相机位置误差影响下，WNCC和NCC算法的陆标匹配结果如图10所示。

图  10  在不同$ {\sigma _r} $影响下，WNCC与NCC的匹配结果

Figure 10.  Landmark matching results of WNCC and NCC under different $ {\sigma _r} $

可以看出，WNCC和NCC的匹配精度均随着陆标位置误差的增大而降低，但WNCC具有更高的精度，其对陆标的位置误差具有更强的稳定性。然而，当${\sigma _r} = 0.1\;{\rm m}$时，WNCC的匹配误差$ [ {\begin{array}{*{20}{c}} {0.06}&{0.04} \end{array}} ]$ pixel略高于NCC$[ {\begin{array}{*{20}{c}} {0.03}&{0.03} \end{array}} ]$ pixel，这是因为较小的系统误差使陆标在图像中的像素变形也较小，所以WNCC的加权过程对精度提升不明显；同时NCC使用了陆标上所有的点（9 801）进行匹配，所以在这种情况下NCC算法具有更高的精度。

同样地，仅考虑相机与陆标之间的相对姿态估计误差，令$ {\sigma _\xi } $由$ 0.075^\circ $开始以$ 0.075^\circ $的间隔均匀增加到$ 1.575^\circ $，$ {\sigma _L} = 0\;{\text{m}} $，$ {\sigma _{L,i}} = 0\;{\text{m}} $且$ {\sigma _r} = 0\;{\text{m}} $。在不同的相机姿态误差影响下，WNCC和NCC算法的陆标匹配结果如图11所示。

图  11  在不同$ {\sigma _\xi } $影响下，WNCC与NCC的匹配结果

Figure 11.  Landmark matching results of WNCC and NCC under different $ {\sigma _\xi } $

可以看出，WNCC与NCC算法的匹配精度均随着姿态误差的增大而降低，但WNCC增长得更加缓慢，体现了其更强的鲁棒性。必须要指出的是，两种算法均对相机的姿态误差较为敏感，当$ {\sigma _\xi } > 2^\circ $时，匹配结果极不稳定；这是因为相机姿态偏差会引起陆标在图像中的像素区域发生严重变形，导致WNCC算法难以选取足够多满足要求的点来执行加权归一化互相关操作。

令WNCC算法选取的陆标点数量由n = 100开始以50的步长均匀增加到n = 1 200，设置系统噪声参数如表1所示。WNCC 利用不同数量的陆标点来进行匹配得到的结果如图12所示。NCC和WNCC在陆标点选取、判断遮挡、图像渲染和图像匹配4个主要步骤中的平均时间消耗，仿真结果如表2所示。

图  12  在不同陆标点数量下，WNCC的匹配结果

Figure 12.  Landmark matching results of WNCC under different number of selected landmark points

由图12可以看出，当n = 100时，陆标匹配结果不可用；但随着陆标点数量的增加，匹配误差迅速降低，当n > 200时，匹配精度几乎没有提升。本算例说明了利用较少的陆标点，通过WNCC算法仍然可以达到较高的匹配精度。需要注意的是，陆标点选取的数量与相机和陆标等参数相关，n > 200仅仅是在本算例条件下获得较高匹配精度的条件。

从表2可以看出，由于进行了陆标点的优化选取，WNCC算法在效率上得到了很大的提升（单核心CPU频率：3.8 GHz）。但在判断遮挡问题上仍有进一步提升的空间。

令相机在陆标坐标系$ {\mathcal{F}^{\text{L}}} $下z方向的坐标由60 m开始以20 m的步长均匀增加到460 m，而x和y方向的坐标始终为0 m，设置系统噪声参数如表1所示。在图像尺度变化的影响下，利用WNCC得到的部分匹配结果如图13所示，其中绿色区域为陆标在图像中的真实位置，红色区域为WNCC估计的位置（下同）。WNCC与NCC的性能比较如图14所示。

图  13  在不同图像尺度下，WNCC的部分匹配结果

Figure 13.  Part of the landmark matching results of WNCC under different image scales

可以看出，两种算法的匹配精度均随着相机与陆标之间距离的降低而降低。因为距离越近，由系统误差引起的陆标在图像中的像素形变越大；但WNCC精度下降的幅度远小于NCC，即对图像尺度变化具有更强的鲁棒性。值得注意的是，由本算例得到的结果并不能表明：相机与陆标的距离越远则匹配精度越高。与此相反，距离过远将导致尺度差异太大，从而导致匹配误差增大甚至算法失效。

图  14  在图像尺度变化影响下，WNCC与NCC的匹配结果

Figure 14.  Landmark matching results of WNCC and NCC under image scale changes

令相机与陆标中心的距离保持为200 m，在${\mathcal{F}^{\rm L}}$的O-xy平面内将相机从$[ {\begin{array}{*{20}{c}} 0&0&{200} \end{array}}]\;{\text{m}} $以4°的间隔移动到$[ {\begin{array}{*{20}{c}} {197.0}&{0.0}&{34.7} \end{array}}]\;{\text{m}} $，即观测视角变化为80°，设置系统噪声参数如表1所示。在相机观测视角变化的影响下，利用WNCC得到的部分匹配结果如图15所示。WNCC与NCC的性能比较如图16所示。

图  15  在不同观测视角下，WNCC的部分匹配结果

Figure 15.  Part of the landmark matching results of WNCC under different viewing geometries

图  16  在观测视角变化影响下，WNCC与NCC的匹配结果

Figure 16.  Landmark matching results of WNCC and NCC under viewing geometry changes

可以看出WNCC算法的匹配精度随着相机与陆标之间视角的增大而降低，而NCC 算法的匹配精度在v方向变化不明显，在u方向逐渐增大；但WNCC算法的误差始终低于NCC 算法，说明其在相机视角变化的情况下具有更高的精度。

令相机在${\mathcal{F}^{\rm L}}$中的位置为$ [ {\begin{array}{*{20}{c}} 0&0&{200} \end{array}} ]\;{\text{m}} $，光照方向在${\mathcal{F}^{\rm G}}$中从$ [ {\begin{array}{*{20}{c}} {0.0}&{ - 0.98}&{ - 0.2} \end{array}}] $以3.1°的步长均匀变化到$ [ {\begin{array}{*{20}{c}} {0.0}&{ - 0.64}&{ - 0.77} \end{array}}] $，即与相机之间的相位角从$ 72.33^\circ $变化到$ 20.26^\circ $，设置系统噪声参数如表1所示。在太阳相位角变化的影响下，利用WNCC得到的部分匹配结果如图17所示。

图  17  在不同相位角下，WNCC的部分匹配结果

Figure 17.  Part of the landmark matching results of WNCC under different phase angles

在不同的太阳相位角下，WNCC与NCC的性能比较如图18所示。可以看出WNCC算法的匹配精度随着相位角的减小而降低，这是因为式（1）描述的光度模型在小相位角的情况下存在较大的误差；而NCC算法的匹配误差在u方向上逐渐增大，在v方向上逐渐降低，这种情况的出现可能与陆标地形相关。但WNCC在光照方向变化的情况下具有更高的精度。

图  18  在相位角变化的影响下，WNCC与NCC的匹配结果

Figure 18.  Landmark matching results of WNCC and NCC under phase angle changes

为应对小天体探测任务当中极端环境对视觉导航可靠性的影响，本文提出了一种陆标鲁棒匹配方法。基于陆标的像素投影模型对其匹配误差进行了详细的理论分析，探讨了相机与陆标之间的相对位姿估计误差对匹配结果的影响；基于误差分析对陆标点进行了优化选取，并设计了WNCC陆标匹配算法；利用高保真度的合成图像序列对比分析了WNCC算法和NCC算法在极端的尺度、视角和光照变化等条件下的匹配性能。数值统计结果表明提出的WNCC算法比NCC算法效率更高、精度更高、鲁棒性更强，更适合应用于小天体探测任务中的视觉导航系统。