李　力，鄧　輝，李　瑧，梅　盈，戴　偉,，楊秋萍,，王　鋒,

(1. 昆明理工大學云南省計算機技術(shù)應(yīng)用重點實驗室，云南 昆明　650500； 2. 中國科學院云南天文臺，云南 昆明　650011；3. 南京大學天文與空間科學學院，江蘇 南京　210000)

光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡(Optical and Near-infrared Solar Eruption Tracer, ONSET)可實現(xiàn)3個通道共4個波段高時間和高空間分辨率的太陽觀測，視場可全日面和局部面切換，是一臺對發(fā)展我國太陽物理觀測研究和空間天氣監(jiān)測具有重大意義的地基太陽望遠鏡[1]。

和所有地基望遠鏡一樣,為克服大氣湍流的影響，獲得接近或者達到衍射極限的圖像[2]，光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡擬采用高分辨圖像重建的方法進行實時數(shù)據(jù)處理，處理過程包括：針對一組原始短曝光圖進行選幀和重建，重建之前進行選幀處理不僅可以提高后續(xù)重建的圖像質(zhì)量，同時在觀測中只保留質(zhì)量較好的圖像，也能減輕原始數(shù)據(jù)的存儲壓力。

顯然，圖像質(zhì)量評價是選幀技術(shù)的關(guān)鍵組成部分。由于在地面觀測中很難得到觀測目標的一幀原始無偏差圖像作為評價參考，因此所有的全參考型和半?yún)⒖夹蛨D像質(zhì)量評價方法如基于結(jié)構(gòu)相似度的方法[3]、均方誤差和峰值信噪比方法[4]均不能用于天文觀測數(shù)據(jù)的選幀處理。而在無參考評價方法中，平均梯度法[5]反映了圖像灰度變化率的大小，并且圖像的細節(jié)反差變化速率衡量的重要標準就是圖像的平均梯度。譜比法[6-7]在太陽高分辨圖像選幀中得到運用，對實驗圖像的要求是圖像要達到足夠多的數(shù)量而且這些圖像在拍攝時間上連續(xù)。

本文的研究工作是針對光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡實時觀測過程中的選幀要求，設(shè)計并實現(xiàn)了一套基于圖形處理器的選幀實時處理模塊，當前模塊實現(xiàn)了平均梯度法和譜比法的高速并行處理。實驗結(jié)果表明，本模塊運行穩(wěn)定可靠，整體實現(xiàn)與性能已經(jīng)可以滿足光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡實時觀測與處理的要求。

1　圖形處理器并行技術(shù)

圖形處理器起初是向中央處理器提供一些基本的操作。但隨著圖形處理器的發(fā)展，已經(jīng)成為具有高并行度的高性能計算處理器。統(tǒng)一計算架構(gòu)(Compute Unified Device Architecture, CUDA)是英偉達公司為圖形處理器增加的一個提高通用計算能力且易用的編程接口。程序員通過使用C語言在統(tǒng)一計算架構(gòu)下編寫程序，就可以在英偉達的圖形處理器下運行*http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html。CUDA工具箱(CUDA Toolkit)包含兩個重要的工具庫，(1)優(yōu)化后的快速傅里葉變換(Fast Fourier Transform)庫CUFFT；(2)線性代數(shù)函數(shù)(Linear Algebra Routine)庫，其中包含了基本線性代數(shù)子程序(Basic Linear Algebra Subprograms, BLAS)CUBLAS。

在中央處理器主機和圖形處理器設(shè)備的數(shù)據(jù)交互中，通常數(shù)據(jù)首先從主機復制到設(shè)備內(nèi)存中，然后由設(shè)備中的線程完成對數(shù)據(jù)預定的操作，并將最后的結(jié)果拷貝返回主機。合理分配圖形處理器設(shè)備中的網(wǎng)格、線程塊、每個塊的線程數(shù)量，采用眾多線程盡可能地達到高并行度，就可以利用圖形處理器強大的計算吞吐量，達到選幀算法的高性能處理的目的。文[8]研究了太陽高分辨圖像重建算法在圖形處理器上的實現(xiàn)。文[6]實現(xiàn)了基于圖形處理器并行計算的Level1級重建程序，使程序在不改變原有處理方式的情況下將重建速度提高近30倍，達到了準實時重建的效果。文[9]在統(tǒng)一計算架構(gòu)下實現(xiàn)了一個Hα全日面云污染實時識別和修復系統(tǒng)。

以上分析可見，圖形處理器并行技術(shù)已經(jīng)在天文領(lǐng)域有諸多應(yīng)用，為開發(fā)光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡觀測數(shù)據(jù)中的實時選幀功能提供了思路，但是針對光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡具體產(chǎn)出的數(shù)據(jù)特點及選幀要求，需要進一步研究開發(fā)能夠適合望遠鏡選幀算法的并行技術(shù)，從而滿足選幀的實時性要求。

2　選幀實時處理模塊的設(shè)計

光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡可實現(xiàn)多通道多波段的觀測，不同的數(shù)據(jù)類型具有不同的圖像特點，不同波段的圖像有可能采用不同的選幀算法。為了滿足觀測過程中的實時選幀要求，有必要設(shè)計并實現(xiàn)一套選幀實時處理模塊。

一般來說，在光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡觀測過程中，數(shù)據(jù)選幀首先是獲取一組待選幀圖像，然后通過圖像質(zhì)量評價方法計算每幀圖像的評價值，依據(jù)評價值按照預期的選幀比例進行挑選，最后將優(yōu)選幀和剔除幀分開存儲，以完成整個選幀過程。

在具體的實現(xiàn)中，如果考慮到中央處理器和圖形處理器的異構(gòu)特點，設(shè)計的選幀算法并行實現(xiàn)模塊流程如圖1。

在實現(xiàn)中針對選幀算法設(shè)計模塊的步驟主要包括以下幾個功能流程：

(1)讀取一組光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡斑點圖(如100幀)到主機內(nèi)存；for parallel implementation

圖1并行實現(xiàn)選幀算法流程圖

Fig.1Flow chart of frame selection algorithm

(2)申請圖形處理器設(shè)備內(nèi)存，拷貝到該組數(shù)據(jù)；

(3)根據(jù)每幀圖像的大小計算設(shè)備執(zhí)行核函數(shù)所需的線程數(shù)；

(4)在圖形處理器中執(zhí)行具體選幀算法邏輯(后面的實驗中具體實現(xiàn)了平均梯度法和譜比法)；

(5)對圖形處理器計算的結(jié)果進行規(guī)約求和，結(jié)果拷貝回主機內(nèi)存；

(6)在主機中，根據(jù)所求每幀評價值，對該組數(shù)據(jù)進行快速排序；

(7)根據(jù)排序結(jié)果將好幀(如評價值前70%)與壞幀分開存放到不同文件目錄。

示意代碼如表1。為今后光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡的發(fā)展，充分考慮到其科學目標的變化、觀測終端的變化，可適應(yīng)待選幀圖像特點，具有較強的開放性和可擴展性。在當前的模塊實現(xiàn)中，在給圖形處理器分配線程數(shù)時，主要依據(jù)待處理每幀圖像的尺寸，使每個像素點可以分配到一個線程，每個線程對應(yīng)處理一個灰度值的計算，整體來看每幀圖像的所有像素點是并行計算和處理的。其中，在圖形處理器中執(zhí)行的選幀算法邏輯部分是選幀的關(guān)鍵，在這個環(huán)節(jié)，可以根據(jù)待選幀圖像的特點和選幀要求，靈活地進行相應(yīng)選幀算法的選取，然后集成到未來的光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡數(shù)據(jù)處理系統(tǒng)中。

表1　選幀實時處理模塊示意代碼Table1　Illustratecodeofframeselectionreal-timeprocessingmodule

總體來看，模塊實現(xiàn)盡可能地優(yōu)化了選幀的執(zhí)行效率，具有較高的并行度。在對每幀選幀評價值規(guī)約求和時，利用圖形處理器中的共享內(nèi)存實現(xiàn)并行規(guī)約思想。首先取得當前圖形處理器中每個線程塊所能分配的最大線程數(shù)，如果所需計算的每幀像素數(shù)小于單個線程塊的最大線程數(shù)，則分配一個線程塊，否則依據(jù)每幀像素數(shù)計算需要分配多少個線程塊。在計算每個線程塊內(nèi)的求和時，由于共享內(nèi)存緩存中的偏移就等于線程索引，因此首先申請共享內(nèi)存緩沖區(qū)所能分配的最大值，然后將線程塊內(nèi)的每個線程計算的值保存在這個緩沖區(qū)，最后對其進行規(guī)約迭代。每個線程塊都按照這樣的策略同時運算，就能高效得到整組選幀數(shù)據(jù)的評價值積分結(jié)果。

3　基于圖形處理器的選幀算法實現(xiàn)

結(jié)合光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡觀測數(shù)據(jù)和選幀算法的特點，目前實現(xiàn)了一種空域方法和一種頻域方法即平均梯度法和譜比法。

3.1　平均梯度

平均梯度法可以體現(xiàn)圖像中物體的邊緣和細節(jié)的強度，求平均梯度值的計算公式如下：

(1)

其中，M和N分別代表圖像的總行數(shù)和總列數(shù)；F(i+ 1,j)為圖像中第i行第j列的灰度值，根號里面代表了F點與其相鄰的下方點和右側(cè)點的灰度值之差的平方和。

從(1)式可以看出，每幀圖像的所有像素點參與計算，必然帶來很大的計算量。在平均梯度法不使用圖形處理器并行加速的情況下，可以通過只計算圖像的部分區(qū)域提升計算速度。因此這里對其進行簡化，即只計算縮小為原圖二分之一的中心區(qū)域。對于一組每幀尺寸為2 560 × 2 160的太陽近全日面圖，選取中心1 280 × 1 080大小的區(qū)域進行計算；對于一組每幀尺寸為1 700 × 1 700的太陽局部面圖，選取中心850 × 850大小的區(qū)域進行計算。

平均梯度法使用圖形處理器并行加速的情況下，針對每幀圖像的全域計算平均梯度值，在算法的并行實現(xiàn)中，每個線程主要參與計算以下幾步，如圖2。關(guān)鍵核函數(shù)的示意代碼如表2。

圖2并行實現(xiàn)平均梯度選幀算法流程圖
Fig.2Flow chart of average gradient frame selection algorithm for parallel implementation

表2　平均梯度選幀算法關(guān)鍵核函數(shù)示意代碼Table2　Illustratecodeofaveragegradientframeselectionalgorithmforkeykernelfunction

(1)判斷當前行列是否越界，計算當前行列對應(yīng)的像素點與其下方點的灰度值之差的平方；

(2)計算當前行列對應(yīng)的像素點與其右側(cè)點的灰度值之差的平方；

(3)將上述兩個平方之和除以2再求均值。

3.2　譜比法

譜比法選幀的主要原理是以目標區(qū)域功率譜的強弱進行選幀評價。主要思路是將每幀圖像的功率譜與該組圖像序列中所有圖像的平均功率譜進行歸一化處理，將相關(guān)目標的量消除，之后在目標圖像的頻段上進行積分，以積分值的大小作為選幀比較的依據(jù)，積分值大的圖像認為質(zhì)量較好。

譜比法需要使用多次傅里葉變換，在不使用圖形處理器并行加速的情況下，只能采用快速傅里葉變換庫(Fastest Fourier Transform in the West, FFTW)計算快速傅里葉變換的操作，由于傅里葉變換較為耗時，對于串行和并行情況下采取截取中心區(qū)域的方法，即對于近全日面圖，確定起始行為540，起始列為640，終止行為1 619，終止列為1 919；對于太陽局部面圖，確定起始行列為425，終止行列為1 274。

采用圖形處理器并行加速的情況下，如圖3，計算譜比選幀評價值的核函數(shù)的處理流程如下：

(1)在截取區(qū)域內(nèi)，每個線程將對應(yīng)的每個點的數(shù)據(jù)由實數(shù)轉(zhuǎn)換為復數(shù)，并且在轉(zhuǎn)換圖像數(shù)據(jù)類型時對每幀大圖進行截取操作；

(2)將截取區(qū)域與窗函數(shù)點乘，即為對截取區(qū)域加窗，同時進行復數(shù)到復數(shù)的傅里葉變換；

(3)計算上一步傅里葉變換后的值，獲取功率譜，同時轉(zhuǎn)換數(shù)據(jù)類型為實數(shù)；

(4)功率譜除以零頻得到歸一化的功率譜，同時將零頻值置為1；

(5)得到該組圖像序列的歸一化功率譜的均值，計算相對功率譜；

(6)相對功率譜與環(huán)函數(shù)點乘并積分，得到在目標頻段上的選幀評價值。在并行算法中采用CUFFT庫代替FFTW庫執(zhí)行快速傅里葉變換的操作。為使每個線程對應(yīng)于圖像中的一個像素點，同時考慮到并行算法中線程數(shù)分配的最優(yōu)化原則，這里根據(jù)截取區(qū)域確定核函數(shù)分配的線程數(shù)。出于優(yōu)化考慮，這里采用CUBLAS庫函數(shù)cublasSscal求平均功率譜，采用cublasSasum計算在目標頻段上的積分。

圖3并行實現(xiàn)譜比法選幀算法流程圖
Fig.3Flow chart of spectral ratio frame selection algorithm for parallel implementation

4　實驗與性能對比

為了驗證本文的算法實現(xiàn)，進行了全面的實驗。測試數(shù)據(jù)共5組，包括1組光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡觀測的Hα太陽近全日面像(每幀尺寸2 560 × 2 160， 占用空間11 MB)和4組具有明顯太陽活動的局部像(前兩組為Hα圖，后兩組為白光圖。每幀尺寸1 700 × 1 700，占用空間6 MB)。在選幀目標上，對每組100幀的近全日面短曝光圖挑選最好的70幀，對每組70幀的局部像挑選50幀。

4.1　平均梯度

在平均梯度法不使用選幀實時處理模塊圖形處理器并行加速的情況下，通過只計算圖像的部分區(qū)域可以提升計算速度。如圖4，在這種提升選幀效率的方式下，部分區(qū)域串行時間明顯低于全圖區(qū)域串行時間，有平均3.7倍的提升。但是，這種方式換來的提升會損失部分有效的圖像質(zhì)量評價信息，影響選幀的準確度。

在平均梯度法使用選幀實時處理模塊圖形處理器并行加速的情況下，針對每幀圖像的全域計算平均梯度值。如圖4，針對全圖計算的并行算法仍然較只計算部分區(qū)域的串行算法速度優(yōu)勢明顯，是全圖區(qū)域串行時間的6～8倍。實驗表明，平均梯度法的選幀實時處理模塊運行穩(wěn)定可靠，實現(xiàn)的并行算法的執(zhí)行效率明顯優(yōu)于原有串行情況下的執(zhí)行效率，有效性顯著。

4.2　譜比法

如圖5，在譜比法使用選幀實時處理模塊圖形處理器并行加速的情況下，對于相同的處理區(qū)域，并行實現(xiàn)譜比算法的執(zhí)行效率明顯優(yōu)于串行情況下的執(zhí)行效率。針對近全日面圖像的選幀總體執(zhí)行時間最快為1.2 s，比原有串行實現(xiàn)提升了7倍；局部面圖最快為0.7 s，平均提升了5倍。實驗表明，集成譜比法后，本文選幀實時處理模塊運行穩(wěn)定可靠，選幀效率提升明顯。

圖4平均梯度法時間比較圖

Fig.4Diagram of average gradient algorithm for time comparison

圖5譜比法時間比較圖

Fig.5Diagram of spectral ratio algorithm for time comparison

5　結(jié)束語

本文提出并實現(xiàn)一套基于圖形處理器的選幀實時處理模塊，對模塊的設(shè)計和實現(xiàn)進行了細致的討論，并已經(jīng)在圖形處理器下實現(xiàn)了平均梯度法和譜比法的高速并行處理。實驗比較得出，模塊取得了較好的加速比，性能已經(jīng)可以滿足光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡實時觀測與處理的要求，為實時高分辨圖像重建打下了堅實的基礎(chǔ)。

隨著光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡數(shù)據(jù)處理和不同科學研究工作的需要，在現(xiàn)有模塊下進一步發(fā)展，研究適合望遠鏡不同類型特點的選幀算法，并實現(xiàn)和集成到未來的光學和近紅外太陽爆發(fā)監(jiān)測望遠鏡數(shù)據(jù)處理系統(tǒng)是下一步研究的內(nèi)容。

參考文獻：

[1]李臻. 太陽低層大氣中小尺度活動的觀測和研究[D]. 南京: 南京大學, 2016.

[2]劉長玉. 選幀技術(shù)在太陽高分辨成像中的應(yīng)用[D]. 昆明: 中國科學院云南天文臺, 2013.

[3]Wang Z, Bovik A C, Sheikh H R, et al. Image quality assessment: from error visibility to structural similarity[J]. IEEE Transaction on Image Processing, 2004, 13(4): 600-612.

[4]鄭艷芳. 太陽高分辨數(shù)據(jù)評價與噴流/日浪的觀測研究[D]. 昆明: 中國科學院云南天文臺, 2014.

[5]趙青, 何建華, 溫鵬. 基于平均梯度和方向?qū)Ρ榷鹊膱D像融合方法[J]. 計算機工程與應(yīng)用, 2012, 48(24): 165-168.

Zhao Qing, He Jianhua, Wen Peng. Image fusion method based on average grads and wavelet contrast[J]. Computer Engineering and Applications, 2012, 48(24): 165-168.

[6]施正, 向永源, 鄧輝, 等. 一米真空太陽望遠鏡Level 1級圖像選幀的GPU實現(xiàn)[J]. 科學通報, 2015, 60(15): 1408-1413.

Shi Zheng, Xiang Yongyuan, Deng Hui, et al. The frame selection for New Vacuum Solar Telescope Level 1 data processing based on GPU[J]. Chinese Science Bulletin, 2015, 60(15): 1408-1413.

[7]劉長玉, 向永源, 金振宇. 太陽色球圖像的選幀處理[J]. 天文研究與技術(shù)——國家天文臺臺刊, 2014, 11(2): 140-144.

Liu Changyu, Xiang Yongyuan, Jin Zhengyu. Image processing with frame selection for Hα solar chromosphere images[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2014, 11(2): 140-144.

[8]向永源. 太陽高分辨高速重建算法的研究[D]. 昆明: 中國科學院云南天文臺, 2016.

[9]張能維, 楊云飛, 李冉陽, 等. Hα全日面云污染實時識別和修復系統(tǒng)[J]. 天文研究與技術(shù), 2016, 13(2): 242-249.

Zhang Nengwei, Yang Yunfei, Li Ranyang, et al. A real-time image processing system for detecting and removing cloud shadows on Hα full-disk solar[J]. Astronomical Research & Technology, 2016, 13(2): 242-249.