基于FPGA的頻域幸運成像算法

黃學明，李彬華,2*，王錦良

(1. 昆明理工大學信息工程與自動化學院，云南 昆明 650500；2. 昆明理工大學云南省計算機應用技術重點實驗室，云南 昆明 650500)

大氣湍流是大氣中一種不規(guī)則的隨機運動，它使望遠鏡的實際分辨率遠低于衍射極限分辨率，是限制地基光學望遠鏡空間分辨率的一個主要因素[1]。地面觀測要獲得接近望遠鏡衍射極限分辨率的圖像，可以使用事后圖像恢復或重建技術，而幸運成像是一種用于去除大氣湍流影響的事后圖像恢復技術[2-3]。幸運成像技術既可以在空域上進行，也可以在頻域上進行，但是空域算法的選圖主要以幀為單位，這樣的選圖方式沒有考慮圖像數據在某些方向的高頻分量信息，在數據處理時，導致圖像數據信息的浪費。因此，文[4]提出了一種基于傅里葉變換的幸運成像算法，該算法主要由圖像預處理、數據選擇和數據疊加3部分組成，并對提出的方法進行了模擬實驗，實驗結果表明，此算法的數據選擇率為10%時，成像效果最佳。文[3]以實測數據進一步證實了頻域幸運成像的有效性，同時指出該算法相對于空域算法在計算量上存在缺陷。文[5]對傳統(tǒng)的頻域幸運成像算法進行改進，在頻域上采取分組數據選擇的方法對傳統(tǒng)的頻域幸運成像算法進行加速，但該算法仍然屬于事后處理的范疇。

上述頻域幸運成像算法都是基于中央處理器，因串行的處理方式，難以實時處理，觀測人員無法在天文觀測過程中準確地了解觀測天體的信息，因此，實時化是幸運成像技術發(fā)展的必然趨勢。近年來，兩大并行的硬件加速器引入了圖像處理領域[6-7]，一個是圖形處理器(Graphics Processing Unit, GPU)，另一個是現場可編程門陣列。雖然兩者都具有強大的并行處理能力，但是現場可編程門陣列是一種半定制型電路，可以根據實際需求編程，而圖形處理器一旦設計完成就固定了，因此，圖形處理器的靈活度遠不如現場可編程門陣列[8]。本文選擇現場可編程門陣列實現頻域幸運成像算法的加速。

本文根據現場可編程門陣列強大的并行性和靈活性，提出了一個基于現場可編程門陣列的頻域幸運成像算法，然后按現場可編程門陣列并行和流水線設計方法，對頻域幸運成像算法用Verilog硬件描述語言(Hardware Description Language, HDL)進行設計，并在Xilinx Spartan-7上進行工程實現和相關測試實驗。

1 適用于現場可編程門陣列的頻域幸運成像算法

傳統(tǒng)的頻域幸運成像算法主要由圖像預處理、頻域內數據選擇和數據疊加3部分組成，在頻域內選擇數據時，需要消耗大量的時鐘資源，導致算法運行速度非常緩慢，屬于事后處理的范疇。本節(jié)提出一種新的頻域幸運成像算法，并在現場可編程門陣列上進行加速，在臺式機上驗證了算法的有效性。

1.1 算法的總體設計

本文提出的基于現場可編程門陣列的頻域幸運成像算法既不同于文[3]和文[4]的傳統(tǒng)頻域幸運成像算法，也不同于文[5]改進的頻域幸運成像算法。文[4]提出在頻域內以幅值為基準對圖像數據進行選擇；文[3]證實了頻域幸運成像有效性的同時，也指出了該算法計算量大的缺陷；文[5]用分組方法對頻域幸運成像算法進行加速，但處理速度仍然非常緩慢。因此，將頻域幸運成像技術在現場可編程門陣列上處理是該算法提高處理速度的有效途徑。

在現場可編程門陣列上進行頻域幸運成像有3個技術難點：

(1)現場可編程門陣列無法對圖像數據直接進行二維快速傅里葉變換。圖像數據在現場可編程門陣列上做二維傅里葉變換與在MATLAB平臺完全不一樣，在MATLAB平臺可以直接對圖像數據做二維變換，而在現場可編程門陣列上只能做一維傅里葉變換。因此，根據二維傅里葉變換的原理，若有一個N行M列的空域圖像數據f(x,y)，二維傅里葉變換后的頻率點數據為F(u,v)，F(u,v)是一個復數，定義為

(1)

對(1)式進行分解，得到

(2)

其中，x，y為空域圖像的行數和列數，取值范圍為x=0, 1, …,M-1；y=0, 1, …,N-1。u，v為頻域圖像的行數和列數，取值范圍為u=0, 1, …,M-1；v=0, 1, …,N-1。由(1)式和(2)式可知，二維傅里葉變換可以采用將二維轉換為兩個一維的算法進行處理，即將二維快速傅里葉變換拆分為行方向的快速傅里葉變換和列方向的快速傅里葉變換，然后按先后次序進行計算。通常，一維快速傅里葉變換算法的行列先后順序對二維快速傅里葉變換運算結果沒有影響，由于本系統(tǒng)是對128 × 128的圖像進行二維快速傅里葉變換，故N，M都等于128。

(2)傳統(tǒng)的頻域幸運成像算法是以幅值為基準選擇數據，即按

(3)

選擇數據。其中，Re[F(u,v)]為頻域數據的實部；Im[F(u,v)]為頻域數據的虛部；M(u,v)為頻域數據的模。

雖然現場可編程門陣列可以方便地進行整數的加減乘除運算，但是在現場可編程門陣列上計算復數幅值的算法復雜且耗時，除非調用Cordic知識產權核(International Property, IP)，但是此核消耗大量的硬件資源，不適合本系統(tǒng)在中小規(guī)模的現場可編程門陣列上實現。因此，本文提出使用頻域數據的實部平方與虛部平方之和代替幅值進行數據選擇，因為實部平方與虛部平方之和與幅值的變化規(guī)律一致，即使用

M2(u,v)=Re[F(u,v)]2+Im[F(u,v)]2

(4)

代替(3)式進行數據選擇。

(3)現場可編程門陣列片的存儲資源有限。由于本系統(tǒng)需要對大量的圖像頻域數據進行存儲，而片上隨機存取存儲器(Random Access Memory, RAM)資源遠遠不足，因此，本系統(tǒng)使用現場可編程門陣列外大容量存儲器進行頻域圖像數據的存儲。

本文提出的適用于現場可編程門陣列的頻域幸運成像算法工作流程為：

(1)圖像輸入，并進行高斯濾波；

(2)濾波后的圖像做行快速傅里葉變換；

(3)行快速傅里葉變換的結果用乒乓操作的方式在隨機存取存儲器中進行轉置；

(4)從隨機存取存儲器中讀出數據做列快速傅里葉變換，并將結果存入第3代雙倍數據率同步動態(tài)隨機存取存儲器(Double-Data-Rate Three Synchronous Dynamic Random Access Memory, DDR3)；

(5)從DDR3中讀出數據，按10%的比例選擇數據；

(6)數據疊加；

(7)對疊加數據做行快速傅里葉逆變換；

(8)行傅里葉逆變換的結果在隨機存取存儲器中進行轉置；

(9)從隨機存取存儲器中讀出數據做列快速傅里葉逆變換；

(10)將列快速傅里葉逆變換的結果分為兩路，一路數據進行二值化，并送入HDMI模塊顯示，另一路數據傳回上位機。

在以上工作流程中，由于現場可編程門陣列具有強大的并行處理能力，第(1)步至第(4)步是并行處理的，第(5)步至第(8)步也是并行處理的。第(1)步圖像數據通過精簡吉比特介質獨立接口(Reduced Gigabit Medium Independent Interface, RGMII)輸入，第(2)步至第(4)步是傅里葉正變換過程，即圖像數據從空域變換到頻域的過程，第(7)步至第(9)步是傅里葉逆變換過程，即圖像數據從頻域變換到空域的過程。

1.2 算法的驗證

為了驗證該頻域幸運成像算法的可行性，同時驗證頻域幸運成像重建的高分辨率圖像比空域幸運成像重建的效果更佳，本文在CPU+MATLAB平臺分別做了兩組實驗，(1)以幅值為基準選擇數據的傳統(tǒng)頻域算法，(2)在MATLAB平臺模擬本文提出的算法在現場可編程門陣列上的實現，并與文[9]的實驗結果進行對比分析。本實驗的硬件平臺是Dell Precision T5500圖像工作站，內存16 GB，中央處理器為Xeon E5620，顯卡為NVIDIA GTX1080，運行的軟件環(huán)境是Windows 7操作系統(tǒng)和MATLAB R2017a。

在MATLAB中對隨機抽取的2 000幀128 × 128像素的圖像進行頻域幸運成像算法處理，根據文[4-5]的研究結果，數據選擇的比例取10%時，頻域幸運成像的效果最好。第1組實驗是對文[4]提出的傳統(tǒng)頻域算法進行復現，其正確性已經過多次驗證，本文不再詳細介紹；第2組實驗是驗證本文提出的頻域幸運成像算法的正確性。這兩組實驗的相關結果如圖1和圖2。

根據上述的實驗結果，直觀上看不出有任何誤差，我們對兩組實驗數據進行分析。如圖1(c)和圖2(c)，由于第2組實驗是模擬算法在現場可編程門陣列上的實現，所以整個運算過程都是對整數進行處理。從圖中數據可以得出，本文提出的適用于現場可編程門陣列的頻域算法跟傳統(tǒng)的算法僅在小數部分存在誤差，即截斷誤差，這是由于現場可編程門陣列只對整數進行處理造成的。同時，為了進一步分析這兩組頻域算法重建高分辨率結果圖像的優(yōu)劣，我們選擇3種客觀的圖像質量評價函數，通過數值計算進行對比分析，其中半峰全寬(Full Width at Half Maxima, FWHM)為天文圖像常用評價指標，數值越小，說明圖像質量越好，其他2種評價函數均為經典的數字圖像清晰度評價指標，數值越大，說明圖像質量越好。表1是3種算法分別對2 000幀圖像處理得到的高分辨率結果圖像的質量評價指標，這3種算法分別是上述兩組頻域幸運成像算法和文[9]空域幸運成像算法。

圖1 傳統(tǒng)算法的結果。(a)三維灰度分布圖；(b)二維灰度圖；(c)峰值數據

圖2 本文提出算法的結果。(a)三維灰度分布圖；(b)二維灰度圖；(c)峰值數據

從表1可以看出，本文提出的頻域算法與傳統(tǒng)的頻域算法相比，半峰全寬是一樣的，騰格拉德(Tenengrad)梯度和拉普拉斯(Laplacia)梯度值略低于傳統(tǒng)的頻域算法，這是現場可編程門陣列運算過程中產生截斷誤差造成的，說明現場可編程門陣列的頻域算法得到的高分辨率圖像質量比傳統(tǒng)算法略差。但是，從表1中本文提出的頻域算法數據和文[9]空域系統(tǒng)的數據對比可知，本文提出的頻域算法得到的圖像質量(清晰度)遠比空域算法得到的圖像質量要好，也說明本文提出的頻域算法是可行有效的。

表1 兩組頻域算法和空域算法的圖像質量評價結果

既然所提算法是可行有效的，那么只要現場可編程門陣列設計合理，就可以縮短處理時間，加速頻域幸運成像，并為實時化提供一種可行的技術方案。所以，在完成上述算法設計后，頻域幸運成像系統(tǒng)的現場可編程門陣列實現就成為下一個關鍵問題。

2 頻域幸運成像算法的現場可編程門陣列設計及系統(tǒng)實現

算法設計及系統(tǒng)的實現通常與所用硬件平臺相關。本文的頻域幸運成像算法的硬件平臺是以Xilinx公司Spartan-7系列的XC7S50芯片為核心的開發(fā)板和一臺中央處理器為Xeon E5620的工作站，設計開發(fā)環(huán)境是Vivado2019.1，采用Vivado平臺的內嵌式邏輯分析儀和ModelSim SE 10.7進行硬件設計的驗證與調試。在系統(tǒng)的研究過程中，將前面提出的算法用Verilog HDL進行了邏輯設計，并在現場可編程門陣列上實現。本文重點介紹系統(tǒng)的總體設計以及FFT_IFFT模塊、數據選擇模塊的算法及現場可編程門陣列實現。

2.1 頻域幸運成像系統(tǒng)的總體設計

本系統(tǒng)主要包含3部分，即數據傳輸、數據運算和數據存儲。由于系統(tǒng)將來要面向實時應用，對數據的傳輸速度要求比較高，因此，系統(tǒng)采用千兆以太網進行數據傳輸。在數據的運算方面，主要涉及二維快速傅里葉變換，而現場可編程門陣列無法對圖像數據直接做二維快速傅里葉變換，故本文使用將二維轉換為兩個一維的方式對天文圖像數據做二維快速傅里葉變換。在數據的存儲方面，系統(tǒng)使用第3代雙倍數據率同步動態(tài)隨機存取存儲器，因為本文是對2 000幀128 × 128像素的圖像進行處理，而現場可編程門陣列的片上隨機存取存儲器資源無法滿足本系統(tǒng)的存儲要求，故采用片外第3代雙倍數據率同步動態(tài)隨機存取存儲器對數據進行儲存，并用Verilog HDL實現系統(tǒng)設計。系統(tǒng)的總體結構框圖如圖3。

圖3 系統(tǒng)的總體結構框圖

由于受到系統(tǒng)使用的現場可編程門陣列硬件平臺片上隨機存取存儲器資源的限制，目前只能使用2 000幀128 × 128像素的短曝光圖像作為原始圖像進行頻域幸運成像處理，并最終以128 × 128像素顯示頻域幸運成像的結果。

2.2 FFT_IFFT算法及實現模塊設計

傳統(tǒng)的空域幸運成像算法是在空域進行選圖和疊加，而頻域幸運成像算法則是在傅里葉圖像的每個頻率點進行數據選擇和疊加，因此，要想實現頻域幸運成像，必須將圖像數據做傅里葉變換，并在頻域進行數據的選擇和疊加。Xilinx公司在Vivado平臺提供了一維快速傅里葉變換的IP核，但由于圖像數據是二維的，因此并不能直接使用一維快速傅里葉變換對空域圖像數據進行運算，因此，本模塊使用將二維轉換為兩個一維的方式對空域圖像做二維傅里葉變換。系統(tǒng)的FFT_IFFT模塊設計結構框圖如圖4。

圖4 2D-FFT模塊的結構框圖

FFT_IFFT模塊的工作流程為(1)將預處理模塊處理后的數據先進行快速傅里葉行變換，再將行變換的結果在RAM1和RAM2中轉置。這里使用兩個深度16 384、寬度48位的簡單雙端口隨機存取存儲器進行乒乓操作(本系統(tǒng)原始圖像數據為16位，經一次傅里葉正變換后，得到的實部和虛部數據各24位，故這里使用寬度48位的隨機存取存儲器進行轉置操作)，即當RAM1寫滿1幀圖像數據時進行讀操作，讀出的數據做快速傅里葉列變換，此時RAM2進行寫操作；當RAM2寫滿1幀圖像數據時進行讀操作，讀出的數據做快速傅里葉列變換，此時RAM1進行寫操作。最后將行列變換后的頻域數據存入第3代雙倍數據率同步動態(tài)隨機存取存儲器進行數據選擇和疊加。(2)將選擇疊加后的數據先進行行IFFT，再將逆變換的結果存入隨機存取存儲器。這里使用一個深度16 384、寬度64位的簡單雙端口隨機存取存儲器進行轉置，即當隨機存取存儲器存完行快速傅里葉逆變換的結果后，再進行列IFFT，即得到疊加后的空域圖像數據，并將此數據分為兩路，一路數據求平均，然后進行重建高分辨率圖像的顯示；另一路經以太網發(fā)送模塊將數據傳回上位機。

2.3 數據選擇算法及實現模塊設計

傳統(tǒng)的數據選擇方式是將所有圖像對應的頻域點的幅值通過排序建立索引關系，選出幅值大小在前10%的頻率點數據，這種數據選擇的方式計算量非常大，并且消耗的硬件資源也很多，即使是Xilinx推出的高端Vertex-7系列現場可編程門陣列也無法滿足這種數據選擇方式所消耗的硬件資源。因此，第一，本文使用選擇數據但不排序的方式進行數據選擇，第二，本文使用頻域數據的實部平方與虛部平方之和代替幅值選擇數據，具體的模塊設計結構框圖如圖5。

圖5 數據選擇模塊的結構框圖

數據選擇的工作流程。首先，從第3代雙倍數據率同步動態(tài)隨機存取存儲器中讀出圖像的頻域數據，讀完第1幀圖像的第1個數據經FIFO(First In First Out)存入隨機存取存儲器，讀出第2幀圖像相對應的數據(即第1個數據)，且以同樣的方式存入隨機存取存儲器，一直到第200幀相對應的數據存入隨機存取存儲器，即從第201幀開始，讀出相對應的數據先存入FIFO，并在該數據存入隨機存取存儲器之前，將隨機存取存儲器中前200相對應點數據的實部平方與虛部平方之和進行比較，找出最小的數值及地址，再將該任意數據的實部平方與虛部平方之和跟前面找出的最小數值進行比較，若大于前面找出的最小值，則將該數據存入前面找出最小值的地址，否則直接舍去。從第202幀到第2 000幀，一直重復第201幀的第1個數據的過程，當第2 000幀相對應的數據存完后，隨機存取存儲器將這200個數據送到疊加模塊進行疊加。其次，將每幀圖像的第2個至第16 384個數據重復前述操作，即完成所有圖像的頻域數據選擇。

3 實驗結果及對比分析

為了驗證本文提出的數據選擇方法在現場可編程門陣列上的實現是正確的，我們使用內嵌式邏輯分析儀對數據選擇模塊的數據進行抓取。為了驗證頻域幸運成像算法在現場可編程門陣列上實現的正確性，本文以同樣的算法在CPU+MATLAB平臺進行模擬實驗，并將此算法在兩個平臺的處理結果進行對比分析。

本文采用2016年10月20日云南天文臺麗江觀測站對天文雙星HDS 70的觀測圖像，共10 000幀，這組圖像的觀測條件和參數見文[10]。對于圖像幀數以及圖像大小的選取，因受現場可編程門陣列硬件資源的限制，我們多次隨機從10 000幀圖像中抽取2 000幀，并裁剪成128 × 128像素的天文目標短曝光圖像。

3.1 頻域幸運成像在現場可編程門陣列上的實驗結果及對比分析

3.1.1 數據選擇算法結果

本系統(tǒng)采用2 000幀圖像每個相對應的頻率點選出200個頻率點的數據選擇方法，即10%的數據選擇比例。為了內嵌邏輯分析儀抓取的數據更加直觀，這里截取前10行的部分數據按10%的比例進行數據選擇，即10行數據選出1行。數據選擇結果如圖6。

從圖6(a)和圖6(b)可以看出，每行相對應數據的實部平方與虛部平方之和最大的那個數被選出，同樣也是幅值最大的那個數被選出，說明內嵌邏輯分析儀抓取的數據結果正確，本文提出的數據選擇算法正確。

圖6 數據選擇的結果。(a)前10行的部分數據；(b)選出的數據

3.1.2 頻域幸運成像在現場可編程門陣列和MATLAB上的結果分析

本系統(tǒng)的硬件平臺是以Xilinx公司Spartan-7系列的XC7S50芯片為核心的開發(fā)板。系統(tǒng)基于現場可編程門陣列實現的頻域幸運成像算法的驗證，采用10%的數據選擇比例。由于在MATLAB和現場可編程門陣列所得的高分辨率圖像不易進行比較，故本系統(tǒng)在現場可編程門陣列和MATLAB上均對算法處理所得的高分辨率圖像做相同的銳化，突出高分辨率圖像中感興趣的目標或區(qū)域，分別得到了如圖7(a)和圖7(b)的銳化圖像。

圖7 頻域幸運成像處理結果。(a)現場可編程門陣列處理結果；(b)MATLAB處理結果

圖7(a)是由以太網回傳到上位機的圖像數據經銳化后調用MATLAB函數顯示的結果，圖7(b)是在MATLAB上進行頻域幸運成像算法處理并做相同銳化后的結果。通過對比可以看出，兩幀圖像完全相同，從而驗證了系統(tǒng)在現場可編程門陣列上實現的頻域幸運成像算法正確。但是，由于現場可編程門陣列只能處理整數，所以在進行快速傅里葉變換時，產生截斷誤差，不過該誤差非常小，通常在一個數以內。為了使數據的對比更加直觀，在重建后的高分辨率圖像數據做平均之前(即200幀圖像疊加后的結果)，我們隨機截取一段現場可編程門陣列處理后的數據(使用內嵌式邏輯分析儀抓取數據)與MATLAB處理后的數據進行對比，發(fā)現大部分數據誤差在200以內，即做完平均后誤差基本在1以內，說明圖像數據在做快速傅里葉變換時產生的截斷誤差對圖像質量的影響很小。隨機截取的數據如圖8(a)和8(b)。

圖8 頻域幸運成像處理結果的部分數據。(a)現場可編程門陣列處理結果；(b)MATLAB處理結果

3.2 時鐘消耗對比

將本文基于現場可編程門陣列頻域幸運成像算法和文[9]提出的基于現場可編程門陣列幸運成像算法(空域)對1 000幀原始圖像進行成像實驗，對比算法運行時間、圖像數據讀取時間和總運行時間，結果如表2。

表2 運行時間

由于頻域幸運成像算法處理的數據量遠比空域幸運成像要多，因此，頻域幸運成像的算法復雜度也比空域幸運成像要高得多，但是從表2可以看出，本文的頻域幸運成像算法在現場可編程門陣列上的處理速度比文[9]的空域系統(tǒng)快0.51 s，說明本算法充分利用了現場可編程門陣列的并行性和靈活性優(yōu)勢。本系統(tǒng)像素讀取速度比文[9]的空域系統(tǒng)快8.4倍，這是因為本系統(tǒng)采用千兆以太網傳輸方式取代SD卡的傳輸方式。本系統(tǒng)的平均運行速度是文[9]空域系統(tǒng)的6.4倍，說明本算法在現場可編程門陣列實現的合理性與正確性。

4 結 論

本文在分析文[4]基于傅里葉變換的幸運成像算法以及文[5]改進的頻域幸運成像算法后，提出了基于現場可編程門陣列的頻域幸運成像算法，并將算法在現場可編程門陣列上實現了加速。首先，本文對算法進行了詳細的概述，并在MATLAB上驗證了該算法是可行有效的，同時驗證了本文的頻域算法在成像效果上優(yōu)于空域算法；其次，本文對算法的總體設計以及二維快速傅里葉變換模塊、數據選擇模塊和數據疊加模塊的設計思路進行了介紹；最后，本文在一塊中小規(guī)模、低成本的現場可編程門陣列開發(fā)板上進行了系統(tǒng)實現，并展示了相應模塊的測試結果與驗證過程。由于MATLAB的代碼效率比較低，不太適合用于天文觀測，而比較合適做算法驗證或實驗結果分析，因此，本文將該頻域系統(tǒng)在現場可編程門陣列上的實驗結果與MATLAB平臺的實驗結果進行了對比，證明了系統(tǒng)在現場可編程門陣列上實現的正確性。雖然本文設計方案還存在不足之處，跟動態(tài)實時化相比還有一定的差距，但是就頻域幸運成像算法在現場可編程門陣列上的具體實現，并在現場可編程門陣列上實現了加速處理，為頻域幸運成像的實時化提供了一種可行的技術方案。

致謝：本文實驗所用雙星HDS 70的短曝光圖像原始數據由中國科學院云南天文臺張西亮提供，特此致謝。