兀 穎，葛 亮，盧曉猛，田健峰，王匯娟，姜曉軍*

(1.中國(guó)科學(xué)院 光學(xué)天文重點(diǎn)實(shí)驗(yàn)室(國(guó)家天文臺(tái))，北京 100101；2.中國(guó)科學(xué)院大學(xué)，北京 100049)

1 引 言

小行星的研究已成為國(guó)際深空探測(cè)領(lǐng)域的研究熱點(diǎn)，其研究方向主要包括：太陽(yáng)系起源和演化；小行星軌道演化機(jī)制；小行星的形成和演化；近地小行星風(fēng)險(xiǎn)防御等[1-2]。小行星的探測(cè)主要分為地面觀測(cè)和空間探測(cè)，利用地基光學(xué)望遠(yuǎn)鏡進(jìn)行小行星觀測(cè)是目前小行星搜尋和性質(zhì)研究的主要手段。其中，卡特琳娜巡天系統(tǒng)(Catalina Sky Survey，CSS) 用于發(fā)現(xiàn)和跟蹤近地天體(Near Earth Objects，NEOs)，該系統(tǒng)中兩臺(tái)巡天望遠(yuǎn)鏡的口徑和視場(chǎng)分別是：1.5 m，15平方度；0.7 m，19.4平方度[3]。ATLAS(Asteroid Terrestrial impact Last Alert System，ATLAS) 主要目的是發(fā)現(xiàn)對(duì)地球有威脅的近地小行星，該系統(tǒng)中目前有兩臺(tái)50 cm 口徑的望遠(yuǎn)鏡，配備有STA1600探測(cè)器，探測(cè)器像元數(shù)為10 560×10 560，視場(chǎng)為29平方度[4-5]。Pan-STARRS(Panoramic Survey Telescope and Rapid Response System，Pan-STARRS) 也開(kāi)展了近地天體的搜尋觀測(cè)，Pan-STARRS1和Pan-STARRS2的口徑和視場(chǎng)都分別為1.8 m， 7平方度[6]。從2014年到2017年，林肯近地小行星研究小組(Lincoln Near-Earth Asteroid Research，LINEAR)使用3.5 m太空監(jiān)視望遠(yuǎn)鏡(Space Surveillance Telescope，簡(jiǎn)稱SST，視場(chǎng)為6平方度)進(jìn)行了大范圍的小行星搜索[7]。

總的來(lái)說(shuō)，小行星地基光學(xué)觀測(cè)設(shè)備向著大視場(chǎng)、大口徑的方向發(fā)展。隨著口徑和視場(chǎng)的增大，探測(cè)能力得到提升，小行星觀測(cè)時(shí)，單幀圖像內(nèi)的恒星數(shù)量也大大增多，需要處理的數(shù)據(jù)量也呈增長(zhǎng)趨勢(shì)。同時(shí)，隨著科學(xué)級(jí)高幀頻CMOS在巡天觀測(cè)中的應(yīng)用，如WFIRST(Weizmann Fast Astronomical Survey Telescope，WFIRST)[8]和TAOS Ⅱ(Trans-Neptunian Automated Occultation Surveys，TAOS)[9]，對(duì)數(shù)據(jù)處理速度的提升有迫切的要求。ATLAS系統(tǒng)中每晚需要處理150 GB原始數(shù)據(jù)，為了實(shí)現(xiàn)數(shù)據(jù)的實(shí)時(shí)處理，該系統(tǒng)中采用兩臺(tái)28核，內(nèi)存為128 GB的機(jī)架式服務(wù)器處理數(shù)據(jù)[4-5]。如何提高數(shù)據(jù)處理速度，是小行星數(shù)據(jù)處理中需要考慮的一個(gè)問(wèn)題。

對(duì)程序進(jìn)行并行化是提升數(shù)據(jù)處理速度的一種有效方式。ATLAS系統(tǒng)采用CPU多線程加快數(shù)據(jù)處理的速度[4-5]。現(xiàn)有的主要并行化處理方法有：基于CPU的并行化、基于GPU的并行化、FPGA硬件并行化。基于CPU的多線程并行計(jì)算，處理速度的提升依賴于CPU核數(shù)，單臺(tái)設(shè)備速率提升空間有限，集群和超級(jí)計(jì)算機(jī)則成本極高； FPGA 屬于純硬件加速，需由專用硬件語(yǔ)言開(kāi)發(fā)，處理速度最快，但算法轉(zhuǎn)換為硬件語(yǔ)言相對(duì)不易，成本較高。

圖像處理器(Graphic Processing Unit，GPU)也稱眾核處理器，不同于CPU，GPU將更多的晶體管用于執(zhí)行單元，在處理單元數(shù)量上遠(yuǎn)超CPU，因此，GPU在處理速度和存儲(chǔ)帶寬上相對(duì)CPU有明顯的優(yōu)勢(shì)。GPU初期只是應(yīng)用于圖像渲染，隨著技術(shù)的發(fā)展和GPU新架構(gòu)的提出，GPU廣泛應(yīng)用于通用計(jì)算領(lǐng)域，近年來(lái)GPU也應(yīng)用于人工智能領(lǐng)域(Artificial Intelligence， AI)。因此，考慮到開(kāi)發(fā)周期、開(kāi)發(fā)成本、速率提升率，本論文采用基于軟件和硬件結(jié)合的GPU加速方式提升小行星光學(xué)觀測(cè)圖像的處理速度。

2 小行星數(shù)據(jù)處理流程

小行星前期數(shù)據(jù)處理的目的是為了獲取其基本信息，用于后期定軌，數(shù)據(jù)處理的流程如圖1所示。

首先對(duì)原始圖像進(jìn)行處理，提取背景；用目標(biāo)圖像減去背景圖像，獲取只包含目標(biāo)的圖像，進(jìn)行目標(biāo)檢測(cè)，提取出目標(biāo)信息；計(jì)算目標(biāo)基本信息；與標(biāo)準(zhǔn)星表進(jìn)行匹配，進(jìn)行恒星證認(rèn)，確認(rèn)靶面坐標(biāo)系和天球坐標(biāo)系之間的轉(zhuǎn)換關(guān)系，計(jì)算轉(zhuǎn)換方程；監(jiān)測(cè)連續(xù)多幀圖像中提取出目標(biāo)的位置變化，檢測(cè)出備選小行星。

圖1 小行星數(shù)據(jù)處理流程

3 基于GPU 的并行化目標(biāo)提取

小行星和恒星的視運(yùn)動(dòng)速度不同[10]，在觀測(cè)圖像中根據(jù)觀測(cè)模式的不同呈現(xiàn)點(diǎn)像或拖長(zhǎng)像。天文圖像處理中一般采用孔徑測(cè)光獲取目標(biāo)的信息，通過(guò)采用不同半徑的圓形孔徑分別圈出目標(biāo)信息(含背景信息)和背景信息，通過(guò)扣除背景信息獲取目標(biāo)信息。但對(duì)于非圓形目標(biāo)，孔徑測(cè)光不再適用。

Source Extractor(SExtractor)是一套從巡天圖像中檢測(cè)天體并提取天體星等、位置等信息的開(kāi)源軟件[11]。SExtractor適用于提取不同形狀和大小的目標(biāo)，且提取算法的魯棒性較好，也適用于提取拖長(zhǎng)星象，可用于小行星光學(xué)觀測(cè)圖像處理。然而，SExtractor算法是串行的，隨著圖像增大和圖像中目標(biāo)數(shù)量的增多，處理速度難以提升。因此本文通過(guò)將SExtractor算法簡(jiǎn)化和并行化后，使其可以在GPU平臺(tái)上運(yùn)行，大幅提高了數(shù)據(jù)處理速度。

3.1 背景提取

采集到的原始圖像包含目標(biāo)信息、天光背景信息、本底等，為了獲取目標(biāo)的準(zhǔn)確信息，需要扣除背景和本底等信息。背景由于構(gòu)成復(fù)雜，在大尺度(102像素量級(jí))和小尺度(10像素量級(jí))上都存在不均勻性，大視場(chǎng)望遠(yuǎn)鏡這一現(xiàn)象尤為明顯，不能用單一值替代。

SExtractor通過(guò)結(jié)合K-σClipping 法和模式估計(jì)法擬合背景[11]，對(duì)圖像進(jìn)行分割，提取子圖特征信息，擬合出背景圖像。整個(gè)計(jì)算過(guò)程計(jì)算量大，且隨著視場(chǎng)增大和探測(cè)能力增強(qiáng)，圖像中包含背景星增多時(shí)，SExtractor算法中原本就耗時(shí)較大的Detection 過(guò)程和 Deblending過(guò)程會(huì)需要更多的處理時(shí)間[6]，因而應(yīng)用串行的SExtractor算法將無(wú)法及時(shí)處理大量高幀頻的數(shù)據(jù)。

3.2 背景提取并行化

分析SExtractor 背景提取算法，有以下特點(diǎn)：背景提取算法中，子圖之間的計(jì)算相關(guān)性小，相互依賴關(guān)系較弱；計(jì)算過(guò)程有輸入大數(shù)據(jù)量，輸出小數(shù)據(jù)量的情況；計(jì)算過(guò)程有輸入小數(shù)據(jù)量，輸出大數(shù)據(jù)量的情況。這些特點(diǎn)均適合應(yīng)用GPU并行化進(jìn)行加速。

對(duì)SExtractor各個(gè)過(guò)程均進(jìn)行并行化，主要并行化方式如下：(1)圖像分割成的N×N子圖，子圖和子圖相互之間沒(méi)有迭代關(guān)系，并行化計(jì)算均值、方差，進(jìn)行高斯統(tǒng)計(jì)；(2)對(duì)濾波算法和樣條插值算法進(jìn)行算法內(nèi)并行化，三次樣條插值的計(jì)算是一個(gè)小數(shù)據(jù)量到大數(shù)據(jù)量的過(guò)程，行和行之間的插值計(jì)算及列和列之間的插值計(jì)算依賴性低，通過(guò)對(duì)計(jì)算過(guò)程進(jìn)行分解，行插值和列差值分別并行計(jì)算。以4 096×4 096圖像為例，圖像分塊示意圖如圖2所示，并行化背景提取處理過(guò)程如圖3所示，背景提取部分計(jì)算量和并行度如表1 所示，算法并行化程度和子圖數(shù)量正相關(guān)。

圖2 圖像分塊示意圖

同時(shí)在統(tǒng)一計(jì)算設(shè)備架構(gòu)CUDA C(Compute Unified Device Architecture，CUDA) 程序編寫(xiě)過(guò)程中采用以下措施提升GPU處理速度：

(1)減少CPU內(nèi)存到GPU內(nèi)存之間的拷貝次數(shù)，僅進(jìn)行一次拷貝；

(2)合理設(shè)置Grid和Block大小，子圖大小設(shè)置為32或者16的倍數(shù)；

(3)在精度保證的前提下，CUDA C中數(shù)據(jù)類型采用float 替代double。

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

表1 4K×4K圖像背景提取部分計(jì)算量和并行度

3.3 目標(biāo)提取并行化

經(jīng)過(guò)背景扣除后的圖像保留了包含目標(biāo)信息的一個(gè)個(gè)獨(dú)立像元，如何將一個(gè)個(gè)獨(dú)立的像元形成目標(biāo)，是目標(biāo)檢測(cè)需要解決的問(wèn)題。

Zhao等[12]對(duì)SExtractor算法實(shí)現(xiàn)了GPU下的并行化，本文中的目標(biāo)提取并行化算法基于該文章中的目標(biāo)提取方法，并進(jìn)行了簡(jiǎn)化。SExtractor的目標(biāo)檢測(cè)算法為了適用于巡天圖像中星系的提取和輪廓重疊的目標(biāo)源的提取，采用了多閾值算法(SExtractor算法中的Deblending過(guò)程和Clearning過(guò)程)進(jìn)行多次迭代。SExtractor算法中的Deblending過(guò)程是為了處理星場(chǎng)中由于視場(chǎng)匹配或自然原因(密近雙星；密集星場(chǎng)，例如星團(tuán)和低銀緯天區(qū)等)帶來(lái)的多目標(biāo)能量分布疊加的問(wèn)題，Clearning過(guò)程是為了剔除Deblending過(guò)程中某些情況下(如初始閾值較低時(shí)具有淺輪廓的目標(biāo)如橢圓星系)檢測(cè)出的虛假目標(biāo)[11]。這兩個(gè)過(guò)程耗時(shí)較多，在整個(gè)算法中時(shí)間占比較大[12]。

對(duì)于小行星圖像，一般較少出現(xiàn)恒星混疊現(xiàn)象，在天文定位時(shí)，若視場(chǎng)內(nèi)參考恒星足夠，則可以不選取疊加恒星進(jìn)行模型擬合，只需從圖像中檢測(cè)出小行星和若干恒星用于天文定位。因此，為了提高處理效率，對(duì)原SExtractor進(jìn)行簡(jiǎn)化。省去Deblending和Clearning過(guò)程；在閾值選取時(shí)，采用方差圖像作為閾值，替代原算法中閾值的多次迭代，僅進(jìn)行一次判斷。簡(jiǎn)化后的目標(biāo)提取算法在整個(gè)算法中的時(shí)間占比降低至34%(文獻(xiàn)[12]中SExtractor串行算法處理4K×4K圖像時(shí)，提取+Deblending過(guò)程+Clearning過(guò)程的時(shí)間占比為82%)。

目標(biāo)提取并行化遵循的原則如下原則：

(1)剔除小于閾值的像素，只記錄滿足閾值要求的數(shù)據(jù)，由于背景像素個(gè)數(shù)遠(yuǎn)遠(yuǎn)大于目標(biāo)像素個(gè)數(shù)，這一措施大大較少了數(shù)據(jù)量；

(2)只對(duì)滿足閾值的像素進(jìn)行索引，進(jìn)行并行化的區(qū)域連通計(jì)算；

(3)根據(jù)標(biāo)記，對(duì)索引進(jìn)行排序，得到連通的段，根據(jù)連通段所占像素個(gè)數(shù)是否滿足要求，剔除不滿足要求的連通段。

簡(jiǎn)化后的目標(biāo)提取算法的計(jì)算量和并行度如表2所示(4K×4K圖像)，最高并行量可達(dá)滿足閾值的像素級(jí)，程序?qū)嶋H運(yùn)行時(shí)的并行度受GPU核數(shù)、內(nèi)存分配及資源調(diào)度等限制。

表2 4K×4K圖像目標(biāo)提取部分計(jì)算量和并行度

本文的目標(biāo)提取算法針對(duì)小行星圖像的特點(diǎn)，對(duì)SExtractor算法經(jīng)過(guò)裁剪和優(yōu)化。為了分析簡(jiǎn)化后精度的變化情況，將本文并行化質(zhì)心提取結(jié)果與CPU下串行SExtractor算法(包含Deblending 和Clean全過(guò)程)的質(zhì)心提取結(jié)果進(jìn)行對(duì)比。如表3所示，本文并行目標(biāo)質(zhì)心提取結(jié)果與串行SExtractor算法質(zhì)心提取結(jié)果相差小于1/10像元，滿足小行星地基光學(xué)觀測(cè)天文定位的精度要求。

表3 質(zhì)心提取誤差

4 匹配算法選取和優(yōu)化

匹配算法通過(guò)對(duì)圖像中提取的恒星與標(biāo)準(zhǔn)星表中的恒星進(jìn)行匹配，證認(rèn)恒星，確認(rèn)靶面坐標(biāo)系和天球坐標(biāo)系之間的轉(zhuǎn)換關(guān)系。

4.1 匹配算法的選取

匹配算法所需處理的數(shù)據(jù)量相對(duì)較小，算法一般迭代較多，因而，匹配算法在CPU下實(shí)現(xiàn)?，F(xiàn)有應(yīng)用較多的匹配軟件有Astrometry.net[13]和Visual Pinpoint[14]等。Astrometry.net離線匹配耗時(shí)較多，Visual Pinpoint 需要專門接口進(jìn)行對(duì)接，這兩款軟件均提供API，但無(wú)法對(duì)內(nèi)部函數(shù)進(jìn)行修改。因此選取開(kāi)源的Match算法[15]。Match算法是基于三角形匹配的匹配算法，算法流程如圖4藍(lán)色部分所示(彩圖見(jiàn)期刊電子版)，在圖像提取的恒星中選取亮星構(gòu)建圖像三角形集合；依據(jù)指向信息和視場(chǎng)大小在星表中搜索恒星，選取亮星構(gòu)建星表三角形集合；在兩個(gè)集合之間依據(jù)最小閾值尋找相似三角形集合；對(duì)相似三角形集合中的頂點(diǎn)進(jìn)行出現(xiàn)次數(shù)的投票；投票率高的點(diǎn)為匹配上的點(diǎn)，利用這些匹配上的點(diǎn)，計(jì)算轉(zhuǎn)換關(guān)系。

4.2 匹配算法的優(yōu)化

在測(cè)試過(guò)程中發(fā)現(xiàn)，在圖像中恒星較多時(shí)，匹配成功率高。但是，當(dāng)視場(chǎng)內(nèi)恒星數(shù)量較少時(shí)，Match算法匹配失敗率較高。表4為一幅匹配失敗的圖像(視場(chǎng)0.18°×0.18°，恒星數(shù)目為17顆(信噪比>5))處理過(guò)程中的投票結(jié)果。經(jīng)過(guò)分析發(fā)現(xiàn)，投票率高的4號(hào)和7號(hào)匹配對(duì)是錯(cuò)誤的匹配對(duì)。Match算法經(jīng)過(guò)投票后只選擇一次備選目標(biāo)，這兩對(duì)錯(cuò)誤的匹配對(duì)會(huì)被選中，參與初始轉(zhuǎn)換關(guān)系的計(jì)算，從而引起計(jì)算結(jié)果出錯(cuò)。

圖4 優(yōu)化后的Match算法流程圖

原Match算法是基于三角形匹配的匹配算法，匹配三角形的選取閾值和三角形邊長(zhǎng)比相關(guān)。視場(chǎng)內(nèi)恒星數(shù)目較少時(shí)，出現(xiàn)大三角形(邊長(zhǎng)較長(zhǎng)的三角形)的概率較高，而大三角形對(duì)閾值不敏感，且原Match算法投票后只進(jìn)行一次備選目標(biāo)選擇，錯(cuò)誤的匹配三角形很容易被選中，從而導(dǎo)致匹配失敗。為了解決稀疏星場(chǎng)匹配失敗的問(wèn)題，嘗試通過(guò)以下兩種方式提高初始轉(zhuǎn)換關(guān)系的準(zhǔn)確性：(1)提高相似三角形選取閾值和投票選取閾值；(2)對(duì)轉(zhuǎn)換關(guān)系的正確性進(jìn)行判斷和迭代。由表4可見(jiàn)錯(cuò)誤的匹配對(duì)投票數(shù)仍然較高，通過(guò)提高投票選取閾值無(wú)法將其排除。且經(jīng)過(guò)測(cè)試發(fā)現(xiàn)，三角形選取閾值提高后，參與投票的星變少，但仍然會(huì)引入錯(cuò)誤匹配對(duì)。

表4 投票結(jié)果

靶面坐標(biāo)系和星表平面坐標(biāo)系之間的關(guān)系符合方程(1)，(x，y)為靶面坐標(biāo)系下圖像中星的坐標(biāo)，(μ，σ) 為平面坐標(biāo)系下星表恒星的坐標(biāo)。Cdelt為探測(cè)器象元比例尺，θ為探測(cè)器旋轉(zhuǎn)角。轉(zhuǎn)換矩陣是一個(gè)酉矩陣：

(1)

5 小行星檢測(cè)

小行星的觀測(cè)根據(jù)觀測(cè)目的的不同主要有兩種觀測(cè)模式：跟蹤模式和凝視模式。采用跟蹤模式時(shí)，小行星呈點(diǎn)像，恒星呈拖長(zhǎng)像；采用凝視模式時(shí)，恒星呈點(diǎn)像，小行星呈拖長(zhǎng)像(見(jiàn)圖5)。無(wú)論是采用那種觀測(cè)模式，小行星和恒星的視運(yùn)動(dòng)速度不同，可以通過(guò)監(jiān)測(cè)連續(xù)多幀圖像(>3)提取出目標(biāo)(恒星和小行星)的位置變化來(lái)檢測(cè)小行星。具體流程如圖6所示。

圖5 凝視模式下小行星在圖像中的移動(dòng)

圖6 小行星檢測(cè)流程圖

6 測(cè)量實(shí)驗(yàn)與結(jié)果

測(cè)量實(shí)驗(yàn)所用圖像為新疆天文臺(tái)南山觀測(cè)站1 m大視場(chǎng)光學(xué)望遠(yuǎn)鏡(簡(jiǎn)稱南山1 m望遠(yuǎn)鏡)采集的小行星圖像，南山1 m望遠(yuǎn)鏡的CCD探測(cè)器像元數(shù)為4 160×4 136，視場(chǎng)為1.3°×1.3°[16]。圖像曝光時(shí)間為180 s，數(shù)據(jù)處理S/N閾值設(shè)置為3，提取目標(biāo)總個(gè)數(shù)為5 000個(gè)。所搭建實(shí)驗(yàn)測(cè)試環(huán)境如下：CPU為i7-6700，計(jì)算機(jī)內(nèi)存大小8 GB，系統(tǒng)為64 bit操作系統(tǒng)；GPU型號(hào)為NVIDIA GeForce GTX 2080Ti，核數(shù)4 352，帶寬616 GB/s，顯存11 GB，開(kāi)發(fā)環(huán)境為CUDA10.0，系統(tǒng)為64 bit操作系統(tǒng)。如表5所示為GPU下和CPU下同一幀圖像100次處理耗時(shí)的平均值，測(cè)試結(jié)果表明基于GPU的并行算法相比于CPU串行算法，速度提升了約17倍。

對(duì)2016年4月27日觀測(cè)的9713小行星進(jìn)行定位，定位結(jié)果和IAU MPC(Minor Planet Center，MPC)發(fā)布的精密星歷對(duì)比(如表6所示)，精度優(yōu)于2″。

表5 數(shù)據(jù)處理算法提速比

表6 9713小行星2016年4月27日觀測(cè)數(shù)據(jù)定位與MPC精密星歷

7 結(jié) 論

為了提升小行星地基光學(xué)觀測(cè)數(shù)據(jù)處理的速度，滿足大口徑大視場(chǎng)小行星地基光學(xué)望遠(yuǎn)鏡數(shù)據(jù)處理實(shí)時(shí)性的要求，本文根據(jù)小行星數(shù)據(jù)處理的特點(diǎn)，對(duì)目標(biāo)提取算法進(jìn)行優(yōu)化和并行化，使其可以在GPU平臺(tái)上實(shí)現(xiàn)軟硬件結(jié)合的加速；對(duì)Match匹配算法進(jìn)行優(yōu)化，提升了匹配算法的準(zhǔn)確性和適用性。基于 NVIDIA GeForce GTX 2080Ti搭建實(shí)驗(yàn)平臺(tái)，處理4K×4K圖像時(shí)間小于200 ms，相比于在CPU(CPU為i7-6700，電腦內(nèi)存大小8 GB，系統(tǒng)為64 bit操作系統(tǒng))下的處理速度提升了約17倍，提高了小行星數(shù)據(jù)處理的效率，實(shí)現(xiàn)了基于GPU的小行星光學(xué)觀測(cè)圖像實(shí)時(shí)處理。目前，4K×4K圖像基于GPU的數(shù)據(jù)處理，對(duì)NVIDIA GeForce GTX 2080Ti顯卡的顯存和核數(shù)的利用仍有冗余，該數(shù)據(jù)處理方法對(duì)于大靶面的圖像，更能發(fā)揮GPU眾核提速的優(yōu)勢(shì)。本方法也適用于其他光學(xué)巡天觀測(cè)圖像處理。但本方法中對(duì)混疊目標(biāo)未進(jìn)行特別處理，在兩個(gè)目標(biāo)混疊嚴(yán)重時(shí)，無(wú)法區(qū)分，未來(lái)將在混疊目標(biāo)的區(qū)分算法上進(jìn)行研究，進(jìn)一步提升算法的適用范圍。

致 謝：本文感謝新疆天文臺(tái)南山1 m望遠(yuǎn)鏡全體工作人員的支持。