兀 颖，葛 亮，卢晓猛，田健峰，王汇娟，姜晓军*

(1.中国科学院 光学天文重点实验室(国家天文台)，北京 100101；2.中国科学院大学，北京 100049)

1 引 言

小行星的研究已成为国际深空探测领域的研究热点，其研究方向主要包括：太阳系起源和演化；小行星轨道演化机制；小行星的形成和演化；近地小行星风险防御等[1-2]。小行星的探测主要分为地面观测和空间探测，利用地基光学望远镜进行小行星观测是目前小行星搜寻和性质研究的主要手段。其中，卡特琳娜巡天系统(Catalina Sky Survey，CSS) 用于发现和跟踪近地天体(Near Earth Objects，NEOs)，该系统中两台巡天望远镜的口径和视场分别是：1.5 m，15平方度；0.7 m，19.4平方度[3]。ATLAS(Asteroid Terrestrial impact Last Alert System，ATLAS) 主要目的是发现对地球有威胁的近地小行星，该系统中目前有两台50 cm 口径的望远镜，配备有STA1600探测器，探测器像元数为10 560×10 560，视场为29平方度[4-5]。Pan-STARRS(Panoramic Survey Telescope and Rapid Response System，Pan-STARRS) 也开展了近地天体的搜寻观测，Pan-STARRS1和Pan-STARRS2的口径和视场都分别为1.8 m， 7平方度[6]。从2014年到2017年，林肯近地小行星研究小组(Lincoln Near-Earth Asteroid Research，LINEAR)使用3.5 m太空监视望远镜(Space Surveillance Telescope，简称SST，视场为6平方度)进行了大范围的小行星搜索[7]。

总的来说，小行星地基光学观测设备向着大视场、大口径的方向发展。随着口径和视场的增大，探测能力得到提升，小行星观测时，单帧图像内的恒星数量也大大增多，需要处理的数据量也呈增长趋势。同时，随着科学级高帧频CMOS在巡天观测中的应用，如WFIRST(Weizmann Fast Astronomical Survey Telescope，WFIRST)[8]和TAOS Ⅱ(Trans-Neptunian Automated Occultation Surveys，TAOS)[9]，对数据处理速度的提升有迫切的要求。ATLAS系统中每晚需要处理150 GB原始数据，为了实现数据的实时处理，该系统中采用两台28核，内存为128 GB的机架式服务器处理数据[4-5]。如何提高数据处理速度，是小行星数据处理中需要考虑的一个问题。

对程序进行并行化是提升数据处理速度的一种有效方式。ATLAS系统采用CPU多线程加快数据处理的速度[4-5]。现有的主要并行化处理方法有：基于CPU的并行化、基于GPU的并行化、FPGA硬件并行化。基于CPU的多线程并行计算，处理速度的提升依赖于CPU核数，单台设备速率提升空间有限，集群和超级计算机则成本极高； FPGA 属于纯硬件加速，需由专用硬件语言开发，处理速度最快，但算法转换为硬件语言相对不易，成本较高。

图像处理器(Graphic Processing Unit，GPU)也称众核处理器，不同于CPU，GPU将更多的晶体管用于执行单元，在处理单元数量上远超CPU，因此，GPU在处理速度和存储带宽上相对CPU有明显的优势。GPU初期只是应用于图像渲染，随着技术的发展和GPU新架构的提出，GPU广泛应用于通用计算领域，近年来GPU也应用于人工智能领域(Artificial Intelligence， AI)。因此，考虑到开发周期、开发成本、速率提升率，本论文采用基于软件和硬件结合的GPU加速方式提升小行星光学观测图像的处理速度。

2 小行星数据处理流程

小行星前期数据处理的目的是为了获取其基本信息，用于后期定轨，数据处理的流程如图1所示。

首先对原始图像进行处理，提取背景；用目标图像减去背景图像，获取只包含目标的图像，进行目标检测，提取出目标信息；计算目标基本信息；与标准星表进行匹配，进行恒星证认，确认靶面坐标系和天球坐标系之间的转换关系，计算转换方程；监测连续多帧图像中提取出目标的位置变化，检测出备选小行星。

图1 小行星数据处理流程

3 基于GPU 的并行化目标提取

小行星和恒星的视运动速度不同[10]，在观测图像中根据观测模式的不同呈现点像或拖长像。天文图像处理中一般采用孔径测光获取目标的信息，通过采用不同半径的圆形孔径分别圈出目标信息(含背景信息)和背景信息，通过扣除背景信息获取目标信息。但对于非圆形目标，孔径测光不再适用。

Source Extractor(SExtractor)是一套从巡天图像中检测天体并提取天体星等、位置等信息的开源软件[11]。SExtractor适用于提取不同形状和大小的目标，且提取算法的鲁棒性较好，也适用于提取拖长星象，可用于小行星光学观测图像处理。然而，SExtractor算法是串行的，随着图像增大和图像中目标数量的增多，处理速度难以提升。因此本文通过将SExtractor算法简化和并行化后，使其可以在GPU平台上运行，大幅提高了数据处理速度。

3.1 背景提取

采集到的原始图像包含目标信息、天光背景信息、本底等，为了获取目标的准确信息，需要扣除背景和本底等信息。背景由于构成复杂，在大尺度(102像素量级)和小尺度(10像素量级)上都存在不均匀性，大视场望远镜这一现象尤为明显，不能用单一值替代。

SExtractor通过结合K-σClipping 法和模式估计法拟合背景[11]，对图像进行分割，提取子图特征信息，拟合出背景图像。整个计算过程计算量大，且随着视场增大和探测能力增强，图像中包含背景星增多时，SExtractor算法中原本就耗时较大的Detection 过程和 Deblending过程会需要更多的处理时间[6]，因而应用串行的SExtractor算法将无法及时处理大量高帧频的数据。

3.2 背景提取并行化

分析SExtractor 背景提取算法，有以下特点：背景提取算法中，子图之间的计算相关性小，相互依赖关系较弱；计算过程有输入大数据量，输出小数据量的情况；计算过程有输入小数据量，输出大数据量的情况。这些特点均适合应用GPU并行化进行加速。

对SExtractor各个过程均进行并行化，主要并行化方式如下：(1)图像分割成的N×N子图，子图和子图相互之间没有迭代关系，并行化计算均值、方差，进行高斯统计；(2)对滤波算法和样条插值算法进行算法内并行化，三次样条插值的计算是一个小数据量到大数据量的过程，行和行之间的插值计算及列和列之间的插值计算依赖性低，通过对计算过程进行分解，行插值和列差值分别并行计算。以4 096×4 096图像为例，图像分块示意图如图2所示，并行化背景提取处理过程如图3所示，背景提取部分计算量和并行度如表1 所示，算法并行化程度和子图数量正相关。

图2 图像分块示意图

同时在统一计算设备架构CUDA C(Compute Unified Device Architecture，CUDA) 程序编写过程中采用以下措施提升GPU处理速度：

(1)减少CPU内存到GPU内存之间的拷贝次数，仅进行一次拷贝；

(2)合理设置Grid和Block大小，子图大小设置为32或者16的倍数；

(3)在精度保证的前提下，CUDA C中数据类型采用float 替代double。

图3 并行化背景提取流程图

表1 4K×4K图像背景提取部分计算量和并行度

3.3 目标提取并行化

经过背景扣除后的图像保留了包含目标信息的一个个独立像元，如何将一个个独立的像元形成目标，是目标检测需要解决的问题。

Zhao等[12]对SExtractor算法实现了GPU下的并行化，本文中的目标提取并行化算法基于该文章中的目标提取方法，并进行了简化。SExtractor的目标检测算法为了适用于巡天图像中星系的提取和轮廓重叠的目标源的提取，采用了多阈值算法(SExtractor算法中的Deblending过程和Clearning过程)进行多次迭代。SExtractor算法中的Deblending过程是为了处理星场中由于视场匹配或自然原因(密近双星；密集星场，例如星团和低银纬天区等)带来的多目标能量分布叠加的问题，Clearning过程是为了剔除Deblending过程中某些情况下(如初始阈值较低时具有浅轮廓的目标如椭圆星系)检测出的虚假目标[11]。这两个过程耗时较多，在整个算法中时间占比较大[12]。

对于小行星图像，一般较少出现恒星混叠现象，在天文定位时，若视场内参考恒星足够，则可以不选取叠加恒星进行模型拟合，只需从图像中检测出小行星和若干恒星用于天文定位。因此，为了提高处理效率，对原SExtractor进行简化。省去Deblending和Clearning过程；在阈值选取时，采用方差图像作为阈值，替代原算法中阈值的多次迭代，仅进行一次判断。简化后的目标提取算法在整个算法中的时间占比降低至34%(文献[12]中SExtractor串行算法处理4K×4K图像时，提取+Deblending过程+Clearning过程的时间占比为82%)。

目标提取并行化遵循的原则如下原则：

(1)剔除小于阈值的像素，只记录满足阈值要求的数据，由于背景像素个数远远大于目标像素个数，这一措施大大较少了数据量；

(2)只对满足阈值的像素进行索引，进行并行化的区域连通计算；

(3)根据标记，对索引进行排序，得到连通的段，根据连通段所占像素个数是否满足要求，剔除不满足要求的连通段。

简化后的目标提取算法的计算量和并行度如表2所示(4K×4K图像)，最高并行量可达满足阈值的像素级，程序实际运行时的并行度受GPU核数、内存分配及资源调度等限制。

表2 4K×4K图像目标提取部分计算量和并行度

本文的目标提取算法针对小行星图像的特点，对SExtractor算法经过裁剪和优化。为了分析简化后精度的变化情况，将本文并行化质心提取结果与CPU下串行SExtractor算法(包含Deblending 和Clean全过程)的质心提取结果进行对比。如表3所示，本文并行目标质心提取结果与串行SExtractor算法质心提取结果相差小于1/10像元，满足小行星地基光学观测天文定位的精度要求。

表3 质心提取误差

4 匹配算法选取和优化

匹配算法通过对图像中提取的恒星与标准星表中的恒星进行匹配，证认恒星，确认靶面坐标系和天球坐标系之间的转换关系。

4.1 匹配算法的选取

匹配算法所需处理的数据量相对较小，算法一般迭代较多，因而，匹配算法在CPU下实现。现有应用较多的匹配软件有Astrometry.net[13]和Visual Pinpoint[14]等。Astrometry.net离线匹配耗时较多，Visual Pinpoint 需要专门接口进行对接，这两款软件均提供API，但无法对内部函数进行修改。因此选取开源的Match算法[15]。Match算法是基于三角形匹配的匹配算法，算法流程如图4蓝色部分所示(彩图见期刊电子版)，在图像提取的恒星中选取亮星构建图像三角形集合；依据指向信息和视场大小在星表中搜索恒星，选取亮星构建星表三角形集合；在两个集合之间依据最小阈值寻找相似三角形集合；对相似三角形集合中的顶点进行出现次数的投票；投票率高的点为匹配上的点，利用这些匹配上的点，计算转换关系。

4.2 匹配算法的优化

在测试过程中发现，在图像中恒星较多时，匹配成功率高。但是，当视场内恒星数量较少时，Match算法匹配失败率较高。表4为一幅匹配失败的图像(视场0.18°×0.18°，恒星数目为17颗(信噪比>5))处理过程中的投票结果。经过分析发现，投票率高的4号和7号匹配对是错误的匹配对。Match算法经过投票后只选择一次备选目标，这两对错误的匹配对会被选中，参与初始转换关系的计算，从而引起计算结果出错。

图4 优化后的Match算法流程图

原Match算法是基于三角形匹配的匹配算法，匹配三角形的选取阈值和三角形边长比相关。视场内恒星数目较少时，出现大三角形(边长较长的三角形)的概率较高，而大三角形对阈值不敏感，且原Match算法投票后只进行一次备选目标选择，错误的匹配三角形很容易被选中，从而导致匹配失败。为了解决稀疏星场匹配失败的问题，尝试通过以下两种方式提高初始转换关系的准确性：(1)提高相似三角形选取阈值和投票选取阈值；(2)对转换关系的正确性进行判断和迭代。由表4可见错误的匹配对投票数仍然较高，通过提高投票选取阈值无法将其排除。且经过测试发现，三角形选取阈值提高后，参与投票的星变少，但仍然会引入错误匹配对。

表4 投票结果

靶面坐标系和星表平面坐标系之间的关系符合方程(1)，(x，y)为靶面坐标系下图像中星的坐标，(μ，σ) 为平面坐标系下星表恒星的坐标。Cdelt为探测器象元比例尺，θ为探测器旋转角。转换矩阵是一个酉矩阵：

(1)

5 小行星检测

小行星的观测根据观测目的的不同主要有两种观测模式：跟踪模式和凝视模式。采用跟踪模式时，小行星呈点像，恒星呈拖长像；采用凝视模式时，恒星呈点像，小行星呈拖长像(见图5)。无论是采用那种观测模式，小行星和恒星的视运动速度不同，可以通过监测连续多帧图像(>3)提取出目标(恒星和小行星)的位置变化来检测小行星。具体流程如图6所示。

图5 凝视模式下小行星在图像中的移动

图6 小行星检测流程图

6 测量实验与结果

测量实验所用图像为新疆天文台南山观测站1 m大视场光学望远镜(简称南山1 m望远镜)采集的小行星图像，南山1 m望远镜的CCD探测器像元数为4 160×4 136，视场为1.3°×1.3°[16]。图像曝光时间为180 s，数据处理S/N阈值设置为3，提取目标总个数为5 000个。所搭建实验测试环境如下：CPU为i7-6700，计算机内存大小8 GB，系统为64 bit操作系统；GPU型号为NVIDIA GeForce GTX 2080Ti，核数4 352，带宽616 GB/s，显存11 GB，开发环境为CUDA10.0，系统为64 bit操作系统。如表5所示为GPU下和CPU下同一帧图像100次处理耗时的平均值，测试结果表明基于GPU的并行算法相比于CPU串行算法，速度提升了约17倍。

对2016年4月27日观测的9713小行星进行定位，定位结果和IAU MPC(Minor Planet Center，MPC)发布的精密星历对比(如表6所示)，精度优于2″。

表5 数据处理算法提速比

表6 9713小行星2016年4月27日观测数据定位与MPC精密星历

7 结 论

为了提升小行星地基光学观测数据处理的速度，满足大口径大视场小行星地基光学望远镜数据处理实时性的要求，本文根据小行星数据处理的特点，对目标提取算法进行优化和并行化，使其可以在GPU平台上实现软硬件结合的加速；对Match匹配算法进行优化，提升了匹配算法的准确性和适用性。基于 NVIDIA GeForce GTX 2080Ti搭建实验平台，处理4K×4K图像时间小于200 ms，相比于在CPU(CPU为i7-6700，电脑内存大小8 GB，系统为64 bit操作系统)下的处理速度提升了约17倍，提高了小行星数据处理的效率，实现了基于GPU的小行星光学观测图像实时处理。目前，4K×4K图像基于GPU的数据处理，对NVIDIA GeForce GTX 2080Ti显卡的显存和核数的利用仍有冗余，该数据处理方法对于大靶面的图像，更能发挥GPU众核提速的优势。本方法也适用于其他光学巡天观测图像处理。但本方法中对混叠目标未进行特别处理，在两个目标混叠严重时，无法区分，未来将在混叠目标的区分算法上进行研究，进一步提升算法的适用范围。

致 谢：本文感谢新疆天文台南山1 m望远镜全体工作人员的支持。