郑雨韵,黄翔宇,毛晓艳
(1.北京控制工程研究所,北京 100094;2.空间智能控制技术全国重点实验室,北京 100094)
太阳系存在着数量庞大的小行星、彗星等小天体,小天体撞击地球虽然概率极低,但一旦发生就会引发巨大灾难,对地球生物的生存造成巨大威胁[1]。针对近地小行星撞击地球的威胁,国内外学者提出了多种消除或降低威胁的策略,其中无核动能高速撞击是目前最为成熟、可行的一种策略。为验证这种动能撞击防御方式的可行性,美国和欧洲航天局(European Space Agency,ESA)开展了双小行星撞击和偏移评估试验(Double Asteroid Redirection Test,DART)[2],中国正在开展相关演示验证任务的论证。为保证探测器能准确地以最大动能撞击到小行星预定区域,需要实现探测器接近、逼近小行星过程的高精度自主导航与控制。
小行星的轨道和轮廓等信息对精确高速撞击至关重要。地面对小行星的测量手段十分有限,小行星星历信息精度较差(误差一般200~1 000 km),需航天器利用导航敏感器对目标小行星成像,根据图像等信息估计小行星的轨道。目标小行星的图像识别是图像导航估计的第一步。小行星尺寸小、亮度弱(一般10星等以上),导航敏感器需要具备暗弱目标成像能力,才能对小行星成像。这导致导航敏感器会同时拍摄到大量暗弱的未知恒星,给目标小行星精确识别带来挑战。传统的目标检测包括基于特征的模板匹配和基于深度学习的目标检测算法。其中基于特征的模板匹配法需要较强的特征以区分目标和背景,如DDIS(Deformable DIversity Similarity)[3];基于深度学习的目标检测在卷积过程中极易丢失细节,卷积后难以保留目标小行星的信息,如Faster R-CNN[4]、YOLOv3[5]。针对恒星抑制,常用星表匹配方法,如三角形匹配算法[6]和基于多特征匹配的快速星图识别算法[7],根据图像上星点与星表中恒星之间的最小角距判断该星点是否为恒星,该方法无法抑制未记录在星表上的恒星,还需要进一步区分未知的背景恒星和目标小行星。
美国已完成多个小行星飞越、着陆、撞击任务,如于2005年发射的“深度撞击号”(Deep Impact)、2006年发射的“新地平线号”(New Horizons)、2016年发射的“欧西里斯–雷克斯号”(OSIRIS-Rex)、2021年发射的DART和“露西号”(Lucy)等。其中Deep Impact光学导航依托AutoNav系统[8]将相机的点扩散函数近似为高斯函数,通过测量已知恒星的点扩散以估计高斯函数的参数,在图像上叠加该高斯函数并用最小二乘法调整位置,以此定位其它恒星和目标物体的位置。KinetX开发的图像处理光学导航软件套件KXIMP为New Horizons和OSIRIS-Rex等任务提供导航服务,根据探测器的相机姿态和视场,将图像中星点像素坐标映射到非线性探测器空间,比对星表从而确认星点是否为恒星[9]。DART的小型机动自主实时导航算法(SMARTNav)根据光斑识别视野中2个最大物体中的较小物体,即撞击目标Didymos-B[10]。Lucy的目标小行星质心定位方法与Stardust类似,采用blobber技术搜索图像中的亮像素,连续亮像素形成光斑,根据光斑大小筛选目标[11]。对于远距离的目标小行星识别,较多使用星表匹配法,该方法还存在无法筛除星表中未记录的恒星、运行时间长等问题。
本文针对目标小行星图像特征弱、未知背景恒星多的特点,研究了利用小行星与背景恒星相对运动特征的目标识别方法。首先采取ORB“改进的具有旋转不变性的FAST和BRIEF算法”(Oriented FAST and Rotated BRIEF,ORB)特征点定位与BEBLID特征点描述结合的方法对帧间图像进行配准,再基于阈值分割法识别星点,逐星点计算相应窗口之间的结构相似性指数,最后完成暗弱目标小行星的检测。
由于探测器在不断运动,探测器的相机坐标系随时间发生变化,帧间序列图像上的恒星背景存在平移、旋转等变换,可使用图像配准方法将前一帧图像的坐标系变换为与当前帧图像相同的坐标系,当前帧图像如图1所示。
配准后的前一帧图像与当前帧图像的同一恒星处于图像的同一位置,而目标小行星与背景恒星存在相对运动,因此可对配准后的前一帧图像与当前帧图像进行逐一像素点比较,从而剔除背景恒星,保留目标小行星。在接近段,探测器与目标小行星的距离较远,目标小行星在图像上所占像素面积极小,目标小行星的相对运动不会影响配准结果。具体配准流程如图2所示。
采用ORB对探测器连续拍摄的目标小行星图像进行特征点定位。ORB是Rublee等[12]在2011年提出的一种计算成本低、计算速度远高于SURF[13]、配准精度与SURF接近的特征提取算法,采用改进的FAST(Features from Accelerated Segments Test)算法[14]进行特征点定位。
1)对每一尺度的图像进行特征点定位。以图像上的像素点p为例,以p为圆心,以3为半径画圆,比较p点灰度值Ip与经过该圆的4个等距像素点的灰度值Ip1、Ip2、Ip3、Ip4。若这2个像素点至少存在一对连续的像素点灰度值与Ip的差值大于阈值h,即该点所处区域灰度值变化较快,该区域通常为某种边缘,则确定p点为图像一关键点,如图3所示。
图3 特征点定位(截取)Fig.3 Feature point localization (part)
2)基于Harris算法对关键点排序。以每个关键点为中心,分别取窗口,计算该窗口移动(u,v)产生的灰度差值,灰度差值越大则该关键点越重要。定义Harris角点响应值为[12]
其中:x、y分别为关键点p的横纵坐标;Ix和Iy分别为该像素点的水平差分和竖直差分;采用高斯函数为窗口内每个像素点赋予权重w;设置k为0.04;选取角点响应值较大的前500个关键点作为用于进一步配准的特征点。
使用ORB特征点定位的当前帧图像如图4所示,只有部分亮度较高的背景恒星被定位,目标行星过于暗弱,不会被视为特征点,因此目标小行星在帧间图像的位移大小不会影响后续配准结果。
图4 使用ORB特征点定位的当前帧图像(截取)Fig.4 Current frame image using ORB for feature point localization (part)
2020年,Suárez等[15]提出了一种有效的学习型二进制描述符BEBLID(Boosted Efficient Binary Local Image Descriptor),是对BELID[16]的一种改进算法,通过比较特征点邻域内一系列像素点对所在方框的平均灰度值大小,得到该关键点的描述符,精度与速度均优于ORB。
BEBLID特征点描述符在特征点邻域内选择用于比较的512对方形区域,中心像素为p1、p2,边长为s。初始化每个弱学习器的阈值T,设置相同的权重,采用Adaboost算法不断迭代p1、p2、s、T,最终输出“1”和“0”,得到长度为512位的二进制描述符。BEBLID比较每一对方形区域的平均灰度值以生成描述符。经过迭代部分方形区域包含了特征点邻域内的星点,即每个特征点的描述符提取的特征不是细微的角点特征,而是该星点邻域内其它星点的亮度与位置关系组成的特征,用于配准的特征为星间结构信息。
其中:γ为学习率,弱学习器hk的权重 αk取相同值。特征提取函数f(x)和阈值T组成弱学习器hk
特征提取函数为特征点邻域内被选择的一对以p1、p2为中心像素,由s×s个像素组成的方形区域的平均灰度值之差
遍历前一帧图像特征点与当前帧图像的特征点,由于BEBLID描述符为二进制描述符,因此计算它们的汉明距离(两个等长字符串之间的汉明距离是两个字符串对应位置的不同字符的个数,即两个二进制描述符之间对应位的值不同的个数)。对于当前帧图像的每一特征点,保留在前一帧图像上与其汉明距离最近的特征点,作为正确匹配的特征点对。
由于探测器的运动,帧间序列图像的恒星背景存在平移和旋转,但恒星之间的相对位置不变,是不同时刻探测器相机坐标系上的同一图像。根据正确匹配的特征点对,对前一帧图像进行仿射变换,可将前一帧图像的恒星转换到当前帧图像上同一恒星的位置上,达到配准的目的。
1)选取匹配最正确的3个特征点对。从所有特征点对中随机选出3个特征点对,根据这3个特征点对计算一个变换矩阵。使用该变换矩阵对所有特征点对进行转换,若转换后的点与原本对应的点相差小于阈值,则视该特征点对为正确匹配点对,并记录正确匹配点对的数量。重复上一过程,直至达到最大迭代次数,保留正确匹配点对数量最多的变换矩阵作为最终的变换矩阵。设置初始最大迭代次数为500,置信度为0.995,阈值为3,在迭代过程中,根据置信度、匹配个数、当前正确匹配点对数量和当前最大迭代次数动态更新最大迭代次数。
2)计算仿射变换矩阵,通用的变换公式为[17]
使用前文选取的匹配最正确的3个特征点对计算6个未知数,得到坐标系之间的转换矩阵。
3)根据仿射变换矩阵对前一帧图像进行变换。将前一帧图像的每一个像素点位置坐标与仿射变换矩阵相乘,获得与当前帧图像同一坐标系的图像,如图5所示。配准后的前一帧图像与当前帧图像上的恒星背景重叠,目标小行星存在位移。
图5 前一帧图像、当前帧图像、配准后的前一帧图像(截取)Fig.5 Previous frame image,current frame image and previous frame image after registration (part)
图像中大部分面积是黑色背景,因此采用先分割再比较的方法,提取当前帧图像的所有星点,以每一星点为中心取包含该星点的小窗口,计算每一窗口截取的当前帧图像与配准后的前一帧图像之间的结构相似性指数,根据结构相似性指数判断该窗口所包含的星点是否为目标星点,完整流程如图6所示。
1)对当前帧图像进行二值化。根据导航敏感器的成像能力和已知的目标小行星星等范围,将图像中像素值大于20的像素点赋予新值255,将其它像素点赋予新值0。
2)对二值化后的当前帧图像进行轮廓检测[18]。轮廓检测的具体步骤:①初始化新边界编号(Newest Border Districts,NBD);②扫描二值化图像,在扫描到相邻像素值突变时,认为该处像素点为外边界起始点,并更新边界编号为NBD+1,将该边界点的像素值赋值为NBD;③检查外边界起始点邻域上的非0像素,若存在非0像素,则将该点的像素值赋值为NBD,若不存在非0像素,则认为该边界点为外边界终止点;④重复第2步的检查方式,检查每个新边界点;⑤重复第2~4步,进行全图扫描。
3)计算轮廓的垂直边界最小矩形。由于恒星与目标小行星的成像均为星点,认为轮廓的垂直边界最小矩形的中心点为星点质心。遍历每个轮廓边界,记录同一边界所有轮廓点的最大、最小横纵坐标,从而获得矩形区域4个顶点的坐标。记录每个星点轮廓矩形的左上角顶点坐标(x,y)、矩形宽度w、矩形高度h,并在图7中绘出。
图7 当前帧图像的所有星点(截取)Fig.7 All star points in the current frame image(part)
结构相似性(Structural Similarity Index)是一种根据亮度、对比度和结构3个指标衡量2幅图像相似度的指标[19]。采用逐星点检测的方法,逐个比对包含各个星点的窗口之间的结构相似性指数。由于在该阶段的图像中目标小行星仅占1~3个像素,w和h为1~2,且计算结构相似性指数的滑动窗口大小至少为3,因此设置包含各星点的窗口以每个星点轮廓矩形的中点(x+w/2,y+h/2)为中心,大小为(w+4)×(h+4)。计算每个窗口截取的当前帧图像与配准后前一帧图像之间的结构相似性指数,记录结构相似性指数最小的窗口,初步认为其包含的星点为目标小行星。
由于恒星配准精度足够高,同一恒星之间的结构相似性指数均大于0.4。目标小行星的帧间位移大小不同,因此存在以下3种情况:①帧间位移远大于配准误差,此时最小的结构相似性指数小于等于0.4,对应的窗口包含目标小行星;②帧间位移小于截取窗口边长,但仍大于配准误差,此时最小的结构相似性指数大于0.4,对应的窗口包含目标小行星;③帧间位移小于配准误差,配准后图像中的目标小行星与当前帧图像中的目标小行星存在重叠,根据截取窗口大小和目标小行星所占像素面积,包含目标小行星的窗口对应的结构相似性指数大于0.85,此时最小的结构相似性指数大于0.4,对应的窗口包含配准误差较大的恒星,而非目标行星。若探测器运动速度较慢,可以选取相隔帧数较多的2帧图像进行检测,即可避免第2、3种情况;若探测器运动速度较快,则相隔帧数较多的2帧图像之间恒星背景重复部分较少,难以实现正确的配准,此时通过以下方法区分。
判断最小的结构相似性指数是否小于阈值0.4,若是,则认为是第1种情况,确定该窗口包含的星点为目标小行星,检测正确;若不是,则进一步判断前一时刻正确检测的目标小行星是否位于当前帧初步检测的目标小行星所在的截取窗口中。若是,则认为是第2种情况,初步检测的目标小行星是正确检测。若不是,则初步检测的目标小行星不是正确检测,使用前一时刻正确检测的目标小行星所在窗口分割配准后的前一帧图像和当前帧图像,计算二者之间的结构相似性指数。若大于阈值0.85,则认为是第3种情况,保留前一时刻正确检测的目标小行星坐标为当前帧图像正确检测的目标小行星坐标。若3种情况均不是,则认为未能检测出目标小行星。
使用上述方法对1 105张图像进行目标小行星检测,正确的检测结果如图8所示。
图8 目标小行星检测结果Fig.8 Target asteroid detection results
用于实验的恒星背景由HIP星表和连续时刻下探测器的姿态生成,其中星等为0~9,视场角为6°以便捕获目标行星。共模拟了1 106张分辨率为1 080×1 080的图像,其中每张图像包含50~100颗恒星和1个半径为1~3个像素的目标小行星。从第2张图像开始,将前一时刻拍摄的图像与当前时刻拍摄的图像进行配准,得到1 105张配准后的图像,从而进行后续的目标小行星检测。
传统的图像配准方法有SURF、AKAZE[20]、ORB等,其中ORB速度最快,AKAZE的精度最高。分别使用传统配准算法和结合ORB特征点定位与BEBLID描述符的算法进行配准,并对使用4种方法配准后的1 105张图像进行目标小行星检测。对检测结果进行评估,若检测的目标小行星与真实的目标小行星像素距离小于等于3,则认为检测正确。对比ORB特征点定位结合BEBLID描述符的方法与传统的图像配准方法SURF、AKAZE、ORB的检测精度和配准时间,使用i5-7200U处理器,2.5 GHz主频,结果如表1所示。
表1 图像配准方法对比Table 1 Comparison of image registration methods
传统的图像配准方法检测正确率极低,难以应用于序列星点图像的配准,因为各个细微星点是占像素面积极小的高斯弥散斑点,难以根据角点特征进行匹配。而ORB与BEBLID结合的方法以星间结构为特征信息进行匹配,在速度和精度上均明显优于传统算法,共检测1 105张星点图像,其中1 103张图像检测正确、2张图像检测错误、0张图像未能检测。
基于图像配准的目标检测常用的方法是进行帧间差分,该方法对配准精度有极大精度要求,否则无法完全抑制恒星,易保留恒星边缘。
另一常用方法是先比较再分割,使用固定大小的滑动窗口对全图进行扫描,计算相应部分之间的结构相似度,得到配准后的前一帧图像与当前帧图像的混叠,最后使用阈值分割法定位目标小行星。该方法存在以下缺点:①图像大部分黑色背景消耗了较长时间;②滑动窗口大小难以选取,较大的滑动窗口难以识别相对运动较小的目标,较小的滑动窗口会消耗较长时间;③在目标行星的相对运动较大时会识别出2个目标,分别为配准后前一帧图像上的目标和当前帧图像上的目标,还需要进一步区分;④在目标行星的相对运动较小时,配准后图像上的目标和当前帧图像的目标在混叠后的图像是部分重叠的,从而无法区分,且识别出的目标坐标存在一定误差。
已发射的小行星探测任务在该阶段多采用星表匹配法。首先根据探测器姿态等信息对星表进行分区。定位图像中的各个星点,根据相机视场角、图像尺寸、相机安装方向、探测器姿态,计算图像上的星点在惯性坐标系下的位置矢量,计算其与星表中各个恒星的夹角。根据视场角和图像大小,设置阈值0.01°,若星表上存在一恒星,其惯性坐标与该星点的惯性坐标之间的夹角小于0.01°,则认为该星点为该恒星。遍历图像上的每一个星点,从而筛除恒星,保留目标小行星。
文中采用先分割再比较的方法与先比较再分割的方法相比,检测速度快,且不需要极高的配准精度即可完成检测。星表匹配法的检测时间受图像中星点数量及相应的星表区域大小影响,检测时间长。由于本文使用星表生成仿真图像,检测结果足够精确,而在实际应用中检测精度会受姿态误差、星表误差等影响,且图像中可能存在星表中未记录的星点,则需要进一步检测目标小行星。3种方法的检测结果对比如表2所示。
表2 目标小行星检测方法对比Table 2 Comparison of target asteroid detection methods
针对小天体撞击任务,采用ORB特征点定位与BEBLID特征点描述结合的方法进行图像配准,相比传统的图像配准方法具有更高的精度和速度;对于与恒星背景存在相对运动的目标小行星,采用了先分割再比较的方法进行检测,相比先比较再分割速度更快,并对配准精度具有更强的鲁棒性。