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

(1.中國科學(xué)院 光學(xué)天文重點實驗室(國家天文臺)，北京 100101；2.中國科學(xué)院大學(xué)，北京 100049)

1 引 言

小行星的研究已成為國際深空探測領(lǐng)域的研究熱點，其研究方向主要包括：太陽系起源和演化；小行星軌道演化機制；小行星的形成和演化；近地小行星風險防御等[1-2]。小行星的探測主要分為地面觀測和空間探測，利用地基光學(xué)望遠鏡進行小行星觀測是目前小行星搜尋和性質(zhì)研究的主要手段。其中，卡特琳娜巡天系統(tǒng)(Catalina Sky Survey，CSS) 用于發(fā)現(xiàn)和跟蹤近地天體(Near Earth Objects，NEOs)，該系統(tǒng)中兩臺巡天望遠鏡的口徑和視場分別是：1.5 m，15平方度；0.7 m，19.4平方度[3]。ATLAS(Asteroid Terrestrial impact Last Alert System，ATLAS) 主要目的是發(fā)現(xiàn)對地球有威脅的近地小行星，該系統(tǒng)中目前有兩臺50 cm 口徑的望遠鏡，配備有STA1600探測器，探測器像元數(shù)為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太空監(jiān)視望遠鏡(Space Surveillance Telescope，簡稱SST，視場為6平方度)進行了大范圍的小行星搜索[7]。

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

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

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

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

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

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

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

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

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

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

3.1 背景提取

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

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

3.2 背景提取并行化

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

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

圖2 圖像分塊示意圖

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

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

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

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

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

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

3.3 目標提取并行化

經(jīng)過背景扣除后的圖像保留了包含目標信息的一個個獨立像元，如何將一個個獨立的像元形成目標，是目標檢測需要解決的問題。

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

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

目標提取并行化遵循的原則如下原則：

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

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

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

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

表2 4K×4K圖像目標提取部分計算量和并行度

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

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

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

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

4.1 匹配算法的選取

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

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

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

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

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

表4 投票結(jié)果

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

(1)

5 小行星檢測

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

圖5 凝視模式下小行星在圖像中的移動

圖6 小行星檢測流程圖

6 測量實驗與結(jié)果

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

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

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

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

7 結(jié) 論

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

致 謝：本文感謝新疆天文臺南山1 m望遠鏡全體工作人員的支持。