胡荣海,黄翔宇,徐 超
(1.北京控制工程研究所,北京 100190;2.空间智能控制技术重点实验室,北京 100190)
小天体探测有助于揭示太阳系的起源与演化,并为行星防御积累技术和经验,具有重大的科学意义和应用价值[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,由中心(特征)点及其周围的陆标点构成。以陆标的中心点为原点建立东-北-天陆标坐标系 FL,如图1所示。利用从不同角度和光照条件下拍摄的小天体图像,SPC技术可以估计出各陆标点在 FL中的高程及其所在表面的法向量和反照率。
图1 由SPC技术生成的陆标及其在小天体中的位置Fig.1 Landmark generated by SPC technology and its position on asteroid surface
陆标匹配过程分为3步:第1步,在给定的观测(相机与陆标的相对位姿)和光照条件下,通过光度模型渲染出陆标在相机中的视觉外观,本文将其定义为渲染图像Iren。第2步,根据观测条件通过投影模型估计陆标上各点在实际拍摄图像I当中的像素坐标,得到映射坐标构成的张量矩阵M。第3步,在I中的一定范围内整体平移M得到新的张量矩阵,并利用插值计算每一个陆标点在I中对应的亮度值,将得到的图像定义为校正图像Irec;重复执行第3步,直到Iren与Irec的NCC值最大,则陆标中心点所对应的像素坐标即为匹配结果。
在真空环境和太阳平行光照射条件下,陆标上一点p通过相机被转换为像素亮度值的过程可以通过光度模型来描述。定义p点所在的表面法向量为n,太阳光入射方向为is、反射方向es(即观测相机的方向),如图2所示。用is表示is与n的夹角(即入射角),es表示es与n的夹角(即反射角),α为es与is的夹角(即太阳相位角),则点p在图像中的像素亮度Ip可建模为
其中:Λ(Is,ti,KC)是太阳光强度Is、相机曝光时间ti和相机参数KC(KC表示将入射光的强度转换为相机CCD像素亮度大小的比例因子)的函数;A为p点的表面反照率;B为图像背景亮度;R(α,is,es)是与图2中所示角度相关的光线反射函数,具体可参考文献[27]。根据式(1)可以计算出所有陆标点的像素亮度值,得到的图像即称为渲染图像Iren。
图2 光度模型中的向量及角度定义Fig.2 Definitions of vectors and angles in photometry
在小天体质心固连坐标系FG下,假设相机的空间位置r∈R3,其姿态可以利用李代数表示为从相机坐标系FC到FG的旋转矩阵可表示为Rc2g=exp(ξ∧),其中ξ∧为ξ 中各元素所构成的反对称矩阵。
陆标的像素投影模型如图3所示。假设陆标的中心点在 FG中的位置为L∈R3,从FL到FG的旋转矩阵为Rl2g,则陆标中心点在FC中的位置V为
图3 陆标的像素投影模型Fig.3 Pixel projection model of landmark
令Li∈R3表示任意陆标点在FL中的位置,其在FC中的坐标Vi为
根据式(3)和式(4)可计算出陆标中心点和各陆标点在像平面内的像素坐标u∈R2,ui∈R2
其中:f为相机的焦距;ux和uy为相机主点在图像坐标系中的像素偏移量。
则L和Li在I中所对应的像素坐标构成了张量矩阵M,通过像素插值可以得到校正图像Irec。
根据陆标匹配过程可知,匹配误差主要来源于光度和像素投影模型误差。然而,光度模型是实际光线反射过程的近似,无法对复杂的多次反射和散射进行准确的建模,故难以对其进行误差分析。本节主要基于陆标的像素投影模型分析陆标匹配误差产生的机制,并为后续陆标点的优化及陆标鲁棒匹配算法的设计提供依据。
考虑到陆标的位置及相机的位姿先验信息均存在噪声,利用上述陆标投影模型将会导致陆标在图像中的投影区域(即M)发生变形,从而造成匹配误差,匹配误差产生的机制如图4所示。由于陆标匹配的目的是在图像I中寻找陆标中心点L所对应的像素坐标,故可以将L在图像中的投影作为参考点,利用陆标点Li与L之间的像素坐标偏移误差来描述陆标在图像中投影区域的变形程度。
图4 陆标匹配误差产生的机制Fig.4 Mechanism of landmark matching error
令上述误差项所引起的di的增量为Δdi,即Li相对于陆标中心点L在图像中像素投影的变形量,则有
式(9)描述了陆标的位置误差和相机的位姿误差如何影响陆标在图像中的变形程度。例如,在完全相同的误差条件下,当相机的位置不同时,Li相对于L的变形程度如图5所示。其中,图5(a)的相机在 FG中的位置为[ -56.7 -89.4 151.1 ]m,图5(b)相机的位置为[ -63.3 -76.8 135.7 ]m。可以看出,距中心点越远,陆标点的变形程度受相机和陆标噪声的影响越大;此外,不同的相机状态也会影响其误差传播。这些发生严重变形的点对应图像I中错误的像素坐标,使得到的校正图像Irec与渲染图象Iren在外观上存在很大的不同,故导致了匹配误差。利用带有误差的参数将陆标投影到图像平面上,得到陆标的实际像素区域与真实像素区域不重合。
图5 陆标上各点相对于中心点的像素变形程度Fig.5 Pixel deformation of landmark
根据以上分析可知,当相机的状态或陆标的位置包含较大噪声时,与陆标中心点相距较远的陆标点在图像中的位置估计结果达到了几个甚至几十个像素误差,利用这些点不但不能提升匹配精度,反而会影响匹配结果;此外,陆标的像素变形情况与相机的状态存在紧密联系,即在相同的噪声情况下,不同相机状态将导致不同的误差分布。因此,动态地选取陆标表面上一些变形程度较小的陆标点生成Iren和Irec并进行匹配将会提高匹配精度,同时根据每个点的变形程度动态地分配不同的权值也能降低陆标变形对匹配结果的影响,从而提升鲁棒性。
假设 ΔL、ΔLi、Δr和Δ ξ的各个分量均为独立同分布的高斯白噪声,其标准差已知,分别为 σL、σL,i、σr和σξ。为了量化描述陆标点的变形程度Δdi在当前相机状态下受 ΔL、ΔLi、Δr和Δ ξ的影响,推导了其协方差并定义了描述其变形的度量因子。根据式(9)可以计算出 Δdi的协方差
其中:I3表示3阶单位矩阵,03表示大小为3×3的零矩阵。定义度量因子δi(单位:pixel)为E{Δdi}的迹的平方根,记作
则 δi表征了在给定的噪声条件下,通过像素投影模型,陆标点Li相对于中心点L在图像中像素的变形程度,其值越大,说明Li越容易受到上述误差项的影响,所以利用Li来进行陆标匹配将会造成更大的匹配误差。
为降低计算复杂度并提高陆标的匹配精度和鲁棒性,设计了陆标点优化选取方法以及WNCC陆标鲁棒匹配方法。
“罗塞塔号”(Rosetta)任务使用的maplet陆标大小为99×99,共由9 801个点组成[25]。理论上利用所有陆标点的信息可以提高匹配结果的精度和稳定性,然而在特定的观测条件及噪声条件下,其中一些点可能存在较大的像素投影误差,不适合用于陆标匹配。为了降低计算复杂度并消除这些点对匹配结果的影响,根据度量因子δ 对陆标点进行优化选取得到集合S为
其中:Tδ(单位:pixel)为度量因子的阈值。
然而,由于地形遮挡,太阳入射光线无法照射的区域将在图像当中呈现出阴影,而其反射光线进入相机视场的陆标点将不可见,如图6所示,其中红色表示在相机中可见的光照区域,黑色表示阴影,灰色表示不可见区域。黑色和灰色区域所覆盖的陆标点都应该从SL中剔除。
图6 地形遮挡造成的阴影与不可见区域Fig.6 Shadows and invisible areas caused by terrain occlusion
为判断陆标点在相机当中是否可见,可以计算该点与相机的连线之间是否与陆标网格模型存在其他交点,例如图6中线段与陆标之间存在交点Lj,故Li不可见,而Lj可见。同理,若陆标点与太阳之间的连线与陆标网格模型存在其他交点,则该点为阴影,如Lp。
此外,若SL中任意两个陆标点在图像中的像素投影距离太近,那么这两点所对应的像素亮度值相近,不具备良好的区分度。可以设置距离阈值Tdis,当两点之间的像素距离小于Tdis时,则剔除其中度量因子较大的一个。
为提高陆标匹配的鲁棒性,根据每一个被选取的点Li∈SL所对应的度量因子δi,定义其在陆标匹配中的权值wi(置信度)
其中:s为尺度因子。
运算符∑表示对向量中的所有元素求和,⊙表示向量中对应元素相乘。在一定范围内平移映射坐标张量M并利用像素插值得到新的Irec,再计算其与Iren之间的归一化互相关值。不断重复该过程,记录其中最大Cw值所对应的像素位置作为陆标在图像I中的匹配结果。
为防止在全图中搜索、匹配陆标,以提高算法效率,需要对搜索范围进行动态估计。根据式(3)推导陆标中心L在图像I中的像素坐标u相对于 ΔL、Δr和Δξ的偏导数,得到雅可比矩阵JL
则u在图像I中像素估计误差的协方差可估计为
定义搜索区域的半径为Rs(pixels)
其中:R0为最小搜索半径,其目的是为了防止Rs过小导致搜索区域无法覆盖陆标在I中的真实位置。
为了在图像中得到亚像素级的陆标定位,通过曲面拟合计算相关值Cw的极值所对应的像素坐标,此处不进行详细讨论。WNCC陆标鲁棒匹配算法的流程如图7所示。
图7 WNCC陆标匹配算法流程Fig.7 Flowchart of WNCC landmark matching algorithm
为验证WNCC 陆标匹配算法的性能,开发了小天体视觉仿真系统(Visual Simulation System,VSS)合成了如图1 所示的陆标(其在FG中的位置为L=[ -38.3 -60.3 101.1 ]Tm,网格大小为99×99,分辨率:0.3 m)在不同距离、不同视角和不同光照等极端条件下的序列图像,并在完全相同的噪声条件下比较了WNCC算法和目前小天体探测任务[7-10]使用的NCC算法的匹配性能。
设计了7组仿真算例,分别讨论了陆标点位置误差、相机与陆标之间的相对位置误差、姿态误差、陆标点数量变化、图像尺度变化、相机视角变化和太阳相位角变化对匹配结果的影响。对于前4组仿真算例,设置光照方向为[ 0 -1 0 ],相机在 FL下的位置为[ 0 0 200 ] m,姿态为[ 180°0°0°](3-1-2旋转序列的欧拉角),且在空间当中保持不变,相机拍摄的图像如图8所示。其中,前3组算例设置WNCC选取n=500个陆标点进行陆标匹配。对于后3组仿真算例,设置标称系统噪声参数如表1所示。每组算例均执行1 000 次蒙特卡洛仿真,得到的统计结果以陆标像素均方根误差(Root Mean Squared Error,RMSE)的形式呈现。值得注意的是,尽管式(1)所描述的光度模型与实际情况存在一定误差,但其对陆标匹配结果的影响非常小,可忽略不计。
图8 相机拍摄的图像(与陆标相距200 m)Fig.8 Image captured by the camera (200 m from the landmark)
表1 标称系统噪声参数Table 1 Nominal noises of system parameters
首先,考虑陆标点Li在FL中的位置误差对匹配结果的影响。令 σL,i从0.05 m开始且以0.05 m的间隔均匀增加到1.05 m,σL=0 m、σr=0 m且σξ=0°。在不同的陆标点位置误差影响下,WNCC和NCC算法的陆标匹配结果如图9所示。
图9 在不同σL,i影响下,WNCC与NCC的匹配结果Fig.9 Landmark matching results of WNCC and NCC under differentσL,i
可以看出,随着 σL,i的增大,WNCC和NCC算法的匹配精度逐渐降低,但WNCC算法具有更高的精度,对陆标点的位置误差具有更强的稳定性。
考虑到陆标中心点L的位置误差与相机位置误差是相对的,故只需要讨论其中一种情况,此处只分析相机位置误差对匹配结果的影响。令 σr由0.1 m以0.2 m的间隔增加到4.1 m,σL=0 m、σL,i=0 m且σξ=0°。在不同的相机位置误差影响下,WNCC和NCC算法的陆标匹配结果如图10所示。
图10 在不同σ r 影响下,WNCC与NCC的匹配结果Fig.10 Landmark matching results of WNCC and NCC under different σr
可以看出,WNCC和NCC的匹配精度均随着陆标位置误差的增大而降低,但WNCC具有更高的精度,其对陆标的位置误差具有更强的稳定性。然而,当σr=0.1 m时,WNCC的匹配误差[ 0.06 0.04 ] pixel略高于NCC[ 0.03 0.03 ] pixel,这是因为较小的系统误差使陆标在图像中的像素变形也较小,所以WNCC的加权过程对精度提升不明显;同时NCC使用了陆标上所有的点(9 801)进行匹配,所以在这种情况下NCC算法具有更高的精度。
同样地,仅考虑相机与陆标之间的相对姿态估计误差,令 σξ由0.075°开始以0.075°的间隔均匀增加到1.575°,σL=0m,σL,i=0m且σr=0m。在不同的相机姿态误差影响下,WNCC和NCC算法的陆标匹配结果如图11所示。
图11 在不同σξ影响下,WNCC与NCC的匹配结果Fig.11 Landmark matching results of WNCC and NCC under different σξ
可以看出,WNCC与NCC算法的匹配精度均随着姿态误差的增大而降低,但WNCC增长得更加缓慢,体现了其更强的鲁棒性。必须要指出的是,两种算法均对相机的姿态误差较为敏感,当σξ>2°时,匹配结果极不稳定;这是因为相机姿态偏差会引起陆标在图像中的像素区域发生严重变形,导致WNCC算法难以选取足够多满足要求的点来执行加权归一化互相关操作。
令WNCC算法选取的陆标点数量由n=100开始以50的步长均匀增加到n=1 200,设置系统噪声参数如表1所示。WNCC 利用不同数量的陆标点来进行匹配得到的结果如图12所示。NCC和WNCC在陆标点选取、判断遮挡、图像渲染和图像匹配4个主要步骤中的平均时间消耗,仿真结果如表2所示。
表2 NCC与WNCC算法不同操作的平均时间消耗Table 2 The average time costs of different operations of NCC and WNCC单位:ms
由图12可以看出,当n=100时,陆标匹配结果不可用;但随着陆标点数量的增加,匹配误差迅速降低,当n> 200时,匹配精度几乎没有提升。本算例说明了利用较少的陆标点,通过WNCC算法仍然可以达到较高的匹配精度。需要注意的是,陆标点选取的数量与相机和陆标等参数相关,n> 200仅仅是在本算例条件下获得较高匹配精度的条件。
图12 在不同陆标点数量下,WNCC的匹配结果Fig.12 Landmark matching results of WNCC under different number of selected landmark points
从表2可以看出,由于进行了陆标点的优化选取,WNCC算法在效率上得到了很大的提升(单核心CPU频率:3.8 GHz)。但在判断遮挡问题上仍有进一步提升的空间。
令相机在陆标坐标系 FL下z方向的坐标由60 m开始以20 m的步长均匀增加到460 m,而x和y方向的坐标始终为0 m,设置系统噪声参数如表1所示。在图像尺度变化的影响下,利用WNCC得到的部分匹配结果如图13所示,其中绿色区域为陆标在图像中的真实位置,红色区域为WNCC估计的位置(下同)。WNCC与NCC的性能比较如图14所示。
图13 在不同图像尺度下,WNCC的部分匹配结果Fig.13 Part of the landmark matching results of WNCC under different image scales
图14 在图像尺度变化影响下,WNCC与NCC的匹配结果Fig.14 Landmark matching results of WNCC and NCC under image scale changes
可以看出,两种算法的匹配精度均随着相机与陆标之间距离的降低而降低。因为距离越近,由系统误差引起的陆标在图像中的像素形变越大;但WNCC精度下降的幅度远小于NCC,即对图像尺度变化具有更强的鲁棒性。值得注意的是,由本算例得到的结果并不能表明:相机与陆标的距离越远则匹配精度越高。与此相反,距离过远将导致尺度差异太大,从而导致匹配误差增大甚至算法失效。
令相机与陆标中心的距离保持为200 m,在 FL的O-xy平面内将相机从[ 0 0 200 ]m以4°的间隔移动到[ 197.0 0.0 34.7 ]m,即观测视角变化为80°,设置系统噪声参数如表1所示。在相机观测视角变化的影响下,利用WNCC得到的部分匹配结果如图15所示。WNCC与NCC的性能比较如图16所示。
图15 在不同观测视角下,WNCC的部分匹配结果Fig.15 Part of the landmark matching results of WNCC under different viewing geometries
图16 在观测视角变化影响下,WNCC与NCC的匹配结果Fig.16 Landmark matching results of WNCC and NCC under viewing geometry changes
可以看出WNCC算法的匹配精度随着相机与陆标之间视角的增大而降低,而NCC 算法的匹配精度在v方向变化不明显,在u方向逐渐增大;但WNCC算法的误差始终低于NCC 算法,说明其在相机视角变化的情况下具有更高的精度。
令相机在 FL中的位置为[ 0 0 200 ]m,光照方向在FG中从[ 0.0 -0.98 -0.2 ]以3.1°的步长均匀变化到[ 0.0 -0.64 -0.77 ],即与相机之间的相位角从72.33°变化到20.26°,设置系统噪声参数如表1所示。在太阳相位角变化的影响下,利用WNCC得到的部分匹配结果如图17所示。
图17 在不同相位角下,WNCC的部分匹配结果Fig.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的匹配结果Fig.18 Landmark matching results of WNCC and NCC under phase angle changes
为应对小天体探测任务当中极端环境对视觉导航可靠性的影响,本文提出了一种陆标鲁棒匹配方法。基于陆标的像素投影模型对其匹配误差进行了详细的理论分析,探讨了相机与陆标之间的相对位姿估计误差对匹配结果的影响;基于误差分析对陆标点进行了优化选取,并设计了WNCC陆标匹配算法;利用高保真度的合成图像序列对比分析了WNCC算法和NCC算法在极端的尺度、视角和光照变化等条件下的匹配性能。数值统计结果表明提出的WNCC算法比NCC算法效率更高、精度更高、鲁棒性更强,更适合应用于小天体探测任务中的视觉导航系统。