基于現場可編程門陣列的幸運成像算法的實現*

趙盼孜，李彬華，毛櫳嘩，陶 勇

(昆明理工大學信息工程與自動化學院，云南 昆明 650051)

地基大口徑光學望遠鏡對天體成像的分辨率受制于大氣湍流[1]。為了消除大氣湍流對目標成像質量的影響，常見且有效的方法是采用自適應光學(Adaptive Optics, AO)技術，以及事后處理方式的圖像復原技術，幸運成像技術就是其中之一。幸運成像技術利用部分短曝光圖像中包含的目標高分辨率的結構信息，在觀測后進行圖像重建，以消除大氣湍流的干擾從而獲得高分辨率圖像[1-2]。二十一世紀初出現的電子倍增電荷耦合器件(Electron Multiplying Charge-Coupled Devices, EMCCD)成像技術具有高分辨率、高讀出速度、低噪聲等優點，適合于曝光時間在毫秒及以上量級的微光成像領域[3]。該技術的出現使幸運成像技術獲得了長足的進步，并在天文觀測領域取得了極大的成功，2013年以前發表的論文已超過200篇[2]。典型的幸運成像系統包括英國劍橋大學的LuckyCam、德國馬普研究所的AstraLux和西班牙加那利群島天文研究所的FastCam，它們在雙星或多星系統中暗伴星的發現、天測天力參數計算等方面，科學產出極為豐富[4-5]。

幸運成像技術是一種簡單可行的圖像復原方法，但由于圖像復原是在觀測完成后的一段時間內進行，其缺點明顯，即天文觀測人員對于所拍攝圖像的實時信息了解不多，難以及時發現并糾正觀測中可能存在的偏差或錯誤。解決此問題的一種辦法是改進算法，增加硬件的處理能力[3]，將幸運成像技術實時或準實時化。

傳統的中央處理器(Central Processing Unit, CPU)對圖像數據的處理采用串行方式，對于大量圖像的幸運重建來說，難以實現實時或準實時化。相較于中央處理器，現場可編程門陣列(Field Programmable Gate Array, FPGA)器件有著先天的并行性和靈活性優勢，能提供強大的并行計算能力，可加速數據或信號的處理，特別是圖像處理[6]，在天文觀測和數據設備中使用越來越多[7]。將現場可編程門陣列用于幸運成像，是將該技術實時或準實時化的一種有效方法。2008年FastCam項目研究過程中，首次將幸運成像技術在現場可編程門陣列上實現，并獲得了實時的高分辨率圖像[1,8]；文[9]將幸運成像算法在現場可編程門陣列和圖形處理器中實現。不過，他們并沒有對算法的現場可編程門陣列實現進行細致的描述。現場可編程門陣列與幸運成像算法結合的研究工作國內尚無報道。

本文對傳統的基于中央處理器和程序設計語言的幸運成像基本算法進行了簡短的分析，結合現場可編程門陣列的特點，提出了一種適合現場可編程門陣列處理的幸運成像方案，用Verilog硬件描述語言(HDL)進行各功能模塊的設計，并將幸運成像算法成功地在現場可編程門陣列上實現，最后利用中國科學院云南天文臺麗江觀測站2.4 m天文望遠鏡拍攝的大量目標短曝光圖像進行硬件算法設計的可行性驗證。

1 基于現場可編程門陣列的幸運成像算法處理流程

幸運成像算法的基本思想是按照一定的評價標準，從一系列短曝光圖像中選出符合評價標準的圖像或區域，之后再對其進行配準疊加處理，以獲得不受或少受大氣湍流影響的高分辨率圖像[2]。本文構建的幸運成像系統在實現過程中按照基本的幸運成像算法的設計要求進行，只是為了適應本系統所用現場可編程門陣列的處理邏輯與資源限制，在不改變幸運成像算法基本處理流程(選圖、配準、疊加)的基礎上將該算法用Verilog HDL重新進行了設計。

1.1 幸運圖像的選取

在一系列短曝光序列圖像中選出 “幸運圖像” 是最終獲得高分辨圖像的關鍵。然而為獲得 “幸運圖像” 必須選取一個合適的像質評價函數。對于空間點源目標圖像來說，像質評價函數通常采用圖像的斯特列爾比(Strehl Ratio, SR)作為評價標準。由于斯特列爾比定義為存在像差時圖像的峰值光強與不存在像差時圖像的峰值光強之比，因此，計算時需要理想的無像差峰值光強。但這個理想的峰值光強在現實中不容易獲得，故通常采用其它的替代方案。一個簡單有效的方法是使用瞬時的斯特列爾比作為圖像的像質評價函數[2,10]，即只對圖像上的峰值光強像素點的灰度值進行斯特列爾比值估算。

對于現場可編程門陣列這類運算資源有限的硬件，計算斯特列爾比只能采用簡單的處理方式。具體來說，在系統實現過程中，直接使用有像差時圖像的峰值光強像素點的灰度值。由于在配準過程中對每幀圖像的峰值光強像素點灰度值的位置有要求，故在系統選圖過程中，以圖像的峰值光強像素點灰度值及其位置共同對圖像進行像質評價。

1.2 圖像的配準與疊加

配準是幸運成像處理中最重要的一環，如果配準不當會直接導致疊加后的圖像模糊不清，分辨率下降，從而無法將暗弱目標圖像顯示出來。由于系統使用恒星的圖像，所以最常用的配準方法有兩種：(1)以整幅圖像中的最大灰度值為中心截取所需的成像區域；(2)以圖像的質心為配準基點。由于受現場可編程門陣列邏輯及資源的限制，系統采用前一種配準方法。

圖像配準完成后，需要對配準圖像進行疊加。在系統中由于所選圖像數量偏小，所以采用直接疊加方法。疊加后所得圖像進行灰度變換(即銳化)可得到易于觀察的高分辨率圖像。

2 幸運成像算法的現場可編程門陣列實現

本文所述的幸運成像算法實現的硬件平臺是以Xilinx公司Spartan-6系列的XC6SLX16芯片為核心開發板。開發環境為ISE Design Suite 14.7，邏輯設計所用HDL是Verilog，采用ChipScope pro 14.7和ModelSim 10.1d進行硬件設計的驗證與調試。

2.1 算法總體設計

由于項目設計的主要目的是在現場可編程門陣列上實現幸運成像算法，因此為讀取圖片像素數據，首先將天文短曝光圖像存儲在開發板外掛的安全數據卡(Secure Digital Card, SD卡)中，以便給算法模塊讀取圖像數據。此外，為滿足算法對速度及片內存儲的要求，同時由于所用芯片邏輯資源和存儲空間的限制，故在數據處理方式上采用數據流，即逐像素處理的方式，但在每個模塊內以及各個模塊間均采用并行的方式對數據進行處理。與此同時在片內數據的存儲上為使處理速度和存儲器之間達到最佳狀態，采用片內的雙端口隨機存儲器(Random Access Memory, RAM)模塊作為數據緩存器。然而，由于所用現場可編程門陣列芯片隨機存儲器資源以及邏輯資源的約束，只能對1 000張128 × 128像素大小的短曝光圖像作選圖處理，并且在配準過程中只能最大截取64 × 64像素大小的圖片進行配準處理，并最終以64 × 64像素大小顯示幸運成像結果。

對于整個算法的硬件構架，為了算法實現過程中修改及調試的方便，采用模塊化的設計方式。主要模塊包括數據的讀寫、保存模塊，短曝光圖像的選圖、配準、疊加模塊以及最終的重建高分辨率圖像顯示模塊，并用Verilog HDL設計實現。基于現場可編程門陣列實現的幸運成像算法的結構框圖如圖1。

圖1 基于現場可編程門陣列的幸運成像算法的結構框圖
Fig.1 Structure block diagram of the lucky imaging algorithm based on FPGA

幸運成像的現場可編程門陣列的工作流程如下：系統上電后，安全數據卡中存儲的天文圖像數據會通過串行外設接口(Serial Peripheral Interface, SPI)以6.25 Mb/s的速度不斷地向安全數據卡頂層模塊傳送，經安全數據卡頂層模塊處理過的數據經先入先出(First-In First-Out, FIFO)緩沖器向第3代雙倍數據率同步動態隨機存取存儲器(Double-Data-Rate Three Synchronous Dynamic Random Access Memory, DDR3)寫數據模塊中寫入(DDR3速度666 Mb/s)，然后通過DDR3寫數據模塊的處理，向現場可編程門陣列芯片外部的DDR3存儲器中寫入圖像數據。在寫數據模塊向DDR3寫入像素數據的同時，選圖模塊不斷地對每幀圖片的像素值進行最大灰度值求解，然后對求出的所有圖片的最大灰度值排序，當排序結束，說明選圖結束，產生選圖結束信號。

其后，向配準模塊發送要截取的原始圖片的序號以及該圖片最大灰度值的位置參數，以便配準模塊向DDR3讀數據模塊發送所選圖片的首像素地址。與此同時，選圖結束信號也會發送給DDR3讀數據模塊，當選圖結束信號和所選圖像首像素地址同時作用于DDR3讀數據模塊時，該模塊將數據從DDR3存儲器中讀出滿足條件的圖片像素進行疊加處理。從圖1看出，配準、DDR3讀出與疊加這3個模塊是串行的，但用HDL設計及現場可編程門陣列實現時，它們是并行的，它們的運行從時間上看是重疊的。

當所選的全部圖片疊加完成后，將結果進行分段線性灰度變換的處理，以增強目標區域的可視性。然后將最終高分辨率圖像的像素數據保存在片內隨機存儲器中。當保存完畢后，通過視頻圖形陣列(Video Graphics Array, VGA)驅動模塊驅動，將灰度圖在視頻圖形陣列顯示器上顯示。

2.2 選圖模塊的設計

在選圖模塊的實現上，系統主要采用最大灰度值求解模塊以及最大灰度值排序模塊共同搭建而成，這兩個子模塊的搭建主要由比較器、計數器以及片內隨機存儲器完成。其中，最大灰度值查找模塊的狀態轉移圖如圖2。

圖2 最大灰度值查找的狀態轉移圖
Fig.2 State transfer diagram for looking up the maximum gray value

在圖2中，myvaild=1表示安全數據卡中的數據可以取出；cnt表示已讀取的一幀圖像的像素數目。其中，cnt=0表示一幀圖像的全部像素已讀完；eve_t=1表示單幀圖像的最大灰度值開始取出；pic_num表示圖像的幀數；wea_max=1表示可以開始向隨機存儲器中寫數據。從2.1節中可知，在本設計中所用的圖像大小為128 × 128，幀數為1 000，所以cnt最大為16 383，pic_num最大為999。

從圖2可以看出，該模塊的主要功能是對從安全數據卡讀出的圖片像素值逐個進行比較，找出每幀圖像的最大灰度值并將其保存在片內隨機存儲器中供排序使用。值得注意的是，在每求出一幀圖像最大灰度值后，在條件觸發時先將其緩存，然后等相應條件到來時將其進行判斷，看是否滿足配準對最大灰度值位置的要求。滿足則在隨機存儲器相應的地址空間中保存該像素值，不滿足則在隨機存儲器相應的地址空間中存入與該最大灰度值相同比特位數的零值。其次，在向隨機存儲器中保存每幀圖像最大灰度值時，將該最大灰度值所在的圖片序號及位置參數一并保存在隨機存儲器相應地址空間，以便進行并行、串行混合排序。

2.3 配準模塊的設計

在配準算法的實現上，采用以最大灰度值為圖片中心、截取64 × 64像素大小進行疊加處理，其設計主要通過選圖模塊中發送的最大灰度值所在圖片序號以及其位置參數，再根據(1)式計算出要截取的圖片首地址，然后發送給DDR3讀數據模塊，用從DDR3存儲器中讀出相應的圖片像素值，供疊加模塊處理，該模塊并不占用片內隨機存儲器資源，只是少量地使用了邏輯片(Slices)資源。

c3_p0_cmd_addr_r=((pic_cnt-1)×65 536)+

(1)

其中，c3_p0_cmd_addr_r表示要讀取的DDR3的地址；pic_cnt表示要讀取的最大灰度值所在的圖片序號；max_cnt表示要讀取的最大灰度值在圖片中的位置參數；INT(max_cnt/27)表示max_cnt除以27之后只取整數部分，在現場可編程門陣列中用移位寄存器實現。

2.4 疊加模塊的設計

對于圖像的疊加主要是通過片內隨機存儲器資源的使用實現的，在疊加模塊設計中占用的片內隨機存儲器個數為12，消耗約37.5%的存儲器資源。該模塊的結構框圖如圖3。

圖3 疊加模塊結構框圖
Fig.3 Block diagram of stacking module

3 實驗結果及分析

為了驗證幸運成像算法在現場可編程門陣列上實現的正確性，擬用MATLAB對相同的幸運成像算法進行編程并將其在計算機上的處理結果與現場可編程門陣列處理結果進行對比。對于現場可編程門陣列上運行幸運成像算法的處理時間問題，只能與計算機的運行時間做比較。所以本文設計了兩個實驗：(1)與MATLAB處理結果的對比實驗；(2)系統自身處理結果與數據的實驗。

幸運成像實驗必須有天文目標短曝光圖像。本文采用2016年10月20日在云南天文臺麗江觀測站對天文雙星HDS 70的觀測圖像，共10 000幀，這組圖像的觀測條件和參數見文[10]。在圖像幀數以及圖像大小的選取上，由于所用的現場可編程門陣列的限制，隨機從10 000幀圖像中抽取1 000幀，并裁剪成128 × 128像素大小的天文目標短曝光圖像進行處理。

3.1 幸運成像算法在CPU + MATLAB上的平臺處理實驗

實驗硬件平臺是Dell Precision T5500圖像工作站、16 GB內存、NVIDIA GTX1050Ti顯卡；運行的軟件環境是Windows 7操作系統、MATLAB R2014a。

將隨機抽取的1 000幀128 × 128像素大小的圖像進行分組，用MATLAB分別對每組圖像進行同樣的算法處理，在實際處理中只截取64 × 64像素大小的目標區域。根據文[10]的研究結果，選圖比取為1%，得到如圖4的高分辨率圖像和三維灰度分布圖。這一結果與文[10]的結果一致。

圖4 1%選圖比下所得的高分辨率圖像。 (a) 二維灰度圖；(b) 三維灰度分布圖
Fig.4 High-resolution image obtained at 1% selection rate. (a) 2D gray image; (b) 3D gray distribution map

3.2 幸運成像算法在現場可編程門陣列上的實驗驗證

硬件平臺是以Xilinx公司Spartan-6系列的XC6SLX16芯片為核心的開發板。基于現場可編程門陣列實現的幸運成像算法設計的驗證，采用1%的選圖比，即從1 000幀選10幀，執行與MATLAB相同幸運成像的算法處理，然后用ChipScope捕捉最終結果圖中最大灰度值周圍部分像素數據值，得到了如圖5(a)圖的經現場可編程門陣列實現的幸運成像算法處理所得的高分辨率圖像部分像素數據值，并與如圖5(b)圖的經MATLAB實現的幸運成像算法處理所得的高分辨率圖像的相同部分像素值做對比。

圖5 分別用現場可編程門陣列與MATLAB實現的幸運成像算法處理所得圖像像素數值

(a) 現場可編程門陣列處理的結果; (b) MATLAB處理的結果
Fig.5 Image pixel values processed by FPGA and MATLAB respectively(a) Results of FPGA processing; (b) Results of MATLAB processing

通過圖5可知，經現場可編程門陣列實現的幸運成像算法處理所得圖像的最終像素數據與在MATLAB上做同樣算法處理所得圖像像素數據完全相同。

從處理效果上看，最終得到的高分辨率圖像像素值較小，直接通過顯示器顯示圖像時，無論在MATLAB上還是現場可編程門陣列做算法處理所得的高分辨率圖像都比較模糊。所以，在現場可編程門陣列和MATLAB上均對算法處理所得的高分辨率圖像做了相同的灰度變換處理——銳化或圖像增強，將最終得到的高分辨率圖像中感興趣的目標或區域突顯出來，分別得到了如圖6(a)、圖6(b)的最終高分辨率圖像。

圖6 基于現場可編程門陣列和MATLAB實現的幸運成像算法處理所得的高分辨率圖像(a) 現場可編程門陣列處理的結果; (b) MATLAB處理的結果
Fig.6 High-resolution images Obtained by lucky imaging algorithm based on the FPGA and the MATLAB(a) Results of FPGA processing; (b) Results of MATLAB processing

圖6(a)是經現場可編程門陣列實現的幸運成像算法處理所得的高分辨率圖像以1 024 × 768的分辨率在擴展圖形陣列(Extended Graphics Array, XGA)顯示器上顯示并通過手機拍攝得到的，在此采用手機拍攝獲得，是因為系統實現所采用的現場可編程門陣列是一個獨立的外部系統，它是脫離個人計算機單獨工作的，所以基于現場可編程門陣列實現的幸運成像算法處理所得的高分辨率圖像要通過單獨的擴展圖形陣列顯示器進行顯示，因此它與個人計算機上MATLAB所得高分辨率結果圖獲取方式不同，它不能像MATLAB一樣直接進行圖像截屏，而只能采用外部的圖像采集設備進行圖像獲取，所以經現場可編程門陣列實現的幸運成像算法處理所得的高分辨率圖像在本文中采用手機拍攝得到。而圖6(b)經MATLAB處理后顯示的最終高分辨率圖像是以1 280 × 768的寬擴展圖形陣列(Wide Extended Graphics Array, WXGA)分辨率顯示的。由于兩圖的顯示器分辨率不同導致縱橫比不同，導致相同圖像在顯示上稍有差異。此外，圖6(a)是手機拍攝獲得的，由于手機拍攝角度及分辨率的不同導致圖像有些許變化。不過，可看出這兩幅圖像是基本相同的，說明經現場可編程門陣列實現的幸運成像算法處理后的結果與經MATLAB實現的幸運成像算法處理后的結果一致，從而說明基于現場可編程門陣列實現的幸運成像算法設計正確。

綜上所述，無論是從最終高分辨率圖像中的像素數據值對比來看，還是從灰度變換后的直觀顯示最終高分辨率圖像對比來看，都說明了基于現場可編程門陣列實現的幸運成像算法設計與實現的正確性。

3.3 現場可編程門陣列實現幸運成像算法的優勢

現場可編程門陣列的并行性和靈活性是其最大的優勢，因此相較于傳統的中央處理器，現場可編程門陣列對同樣算法的處理速度肯定快。但對于研究的幸運成像算法來說，具體快多少，必須通過實驗說明。為此，分別在現場可編程門陣列和CPU + MATLAB上對相同幸運算法的運算速度進行了測試。實驗過程中，統計1 000幀圖像的處理時間，選圖比采用1%。現場可編程門陣列及CPU + MATLAB平臺上均采用2.1節中介紹的算法實現。由此得到了不同平臺下相同算法的平均運行時間：對于現場可編程門陣列是2.45 s，對于CPU + MATLAB是46.93 s；運行所得結果(高分辨率圖像)分別如圖6(a)、圖6(b)。顯然，相同的算法在不同的平臺處理下，得到相同的高分辨率圖像，現場可編程門陣列上算法處理速度比CPU + MATLAB平臺上算法處理速度快約19倍。如果采用更先進的現場可編程門陣列，速度還將更快。

不過受制于安全數據卡的速度，總體運行時間并沒有明顯的優勢。目前正在設計基于USB3和千兆以太網的高速數據接口，以加速整個幸運成像處理的速度。另外，采用處理能力更強、邏輯資源更多的現場可編程門陣列，可以更多地引入并行計算、減少串行處理過程，充分發揮現場可編程門陣列速度的優勢。

4 結 論

本文首先簡述了幸運算法的基本思想及處理流程，根據所用現場可編程門陣列的硬件資源，用Verilog HDL設計了幸運成像算法的選圖、配準、疊加等關鍵模塊，詳細介紹了這些模塊的實現方案，并在現場可編程門陣列上實現了這一幸運成像系統；然后，通過將現場可編程門陣列處理的結果與個人計算機上用MATLAB處理得到的結果進行對比，證明了此系統設計和實現過程的正確性；最后，通過對現場可編程門陣列與CPU + MATLAB平臺下實現相同算法所需的處理時間進行對比分析，說明了現場可編程門陣列在處理速度上快19倍的優勢。雖然本設計方案還有一些可以改進之處，但就幸運成像算法在現場可編程門陣列上的具體實現，為構建實時化的高速幸運成像系統探索了一種可行的技術方法。

致謝：感謝中國科學院云南天文臺張西亮、季凱帆、金振宇為本文工作提供了原始天文圖像。