楊秋萍，宋振旗，鄧　輝，王　鋒

(昆明理工大學云南省計算機技術應用重點實驗室，云南 昆明　650500)

地基望遠鏡成像受大氣湍流的影響導致圖像的閃爍和抖動從而影響圖像分辨率[1]，而高分辨率太陽觀測圖像是研究太陽物理的必備條件，這就要求望遠鏡能夠實現(xiàn)衍射極限分辨率的成像[2]。

目前事后圖像處理主要包括盲退卷積、相位差法和斑點統(tǒng)計重建等方法*http://adsabs.harvard.edu/abs/1980itek.rept.....H。其中斑點統(tǒng)計重建方法的算法本質是基于一組含有高頻信息的短曝光圖像的統(tǒng)計計算重建目標的高頻信息[3]，如常用的Knox-Thompson(K-T)算法。但是由于觀測數(shù)據(jù)量大、計算復雜等因素的制約，事后圖像重建過程耗時嚴重，無法達到實時觀測的需求。目前國內外先進的望遠鏡通常使用高性能計算機或者大規(guī)模的計算機集群，通過并行計算進行高分辨圖像重建，比如美國大熊湖天文臺[4]采用80核的計算集群。但是大規(guī)模計算集群的成本太高，集群規(guī)模小處理速度滿足不了實時處理的需求，因此迫切需求一種更高效的方法提高計算效率。

近年來圖形處理器通用計算技術發(fā)展迅猛，在許多復雜的計算任務中得到廣泛應用。本文對CPU-GPU混合架構下的高分辨圖像重建進行研究，系統(tǒng)研究了使用MPI-CUDA混合并行技術實現(xiàn)基于K-T算法的高分辨圖像重建的方法，最終實現(xiàn)了一套在單機環(huán)境下的高分辨太陽圖像重建系統(tǒng)。

1　混合編程模型

1.1　CUDA簡介

圖形處理器最早是為提高計算機圖形圖像顯示的實時性和高效性而誕生的，目前它已發(fā)展成為了高密集度、高并行度、多線程、擁有強大計算能力和極高存儲器帶寬的多核處理器[5]，其計算能力遠遠超過了中央處理器，如NVIDIATeslaK80的單精度計算能力達到了7 TFLOPS。

CUDA(Compute Unified Device Architecture)是一種通用的圖形處理器并行計算架構，該架構利用圖形處理器的處理能力能夠大幅提升計算性能，開發(fā)人員可以直接使用C語言等在該架構下編寫程序。運行在圖形處理器上的CUDA并行計算函數(shù)稱為kernel(內核函數(shù))，一個完整的CUDA程序由一系列的kernel函數(shù)和主機端的串行處理部分共同組成。一個kernel函數(shù)以線程網(wǎng)格的形式組織，每個網(wǎng)格由若干個線程塊組成，每個線程塊又由若干個線程組成。一個kernel以線程塊為單位執(zhí)行，只有在同一個線程塊中的線程才可以互相通信。CUDA為一些復雜的計算提供了很多應用程序接口，例如本實驗用到的CUBLAS和CUFFT庫，CUBLAS庫提供了大量對矩陣的基礎性操作，CUFFT庫提供了一到三維各種數(shù)據(jù)類型的傅里葉變換函數(shù)。

1.2　MPI-CUDA混合編程

MPI是一種消息傳遞協(xié)議，它通過消息傳遞實現(xiàn)不同進程之間的通信，從而達到并行計算的目的。MPI在許多復雜的計算任務中得到廣泛應用，明安圖射電頻譜日象儀就使用MPI實現(xiàn)了高性能的UVFITS數(shù)據(jù)合成[6]。但是MPI的并行僅僅是計算任務的并行，它將計算任務分成多個子任務同時執(zhí)行，從而減小整體的計算時間。如果子任務是數(shù)據(jù)密集型的計算，中央處理器的計算能力依然不能滿足實時計算的需求。使用MPI實現(xiàn)對任務的劃分，使用CUDA實現(xiàn)對子任務的計算加速，這就構成了MPI-CUDA混合并行的模式。在中央處理器和圖形處理器架構下，這種模式可以結合兩者的優(yōu)勢，提高科學計算的效率。

2　斑點統(tǒng)計重建方法

斑點統(tǒng)計重建方法是一種以短曝光圖像的統(tǒng)計信息為基礎的方法，這種方法通過對目標斑點圖像進行某種形式的統(tǒng)計處理得到高分辨重建像[7]。按重建算法主要分為頻域重建和空域重建兩類。頻域重建的代表方法有基于二階統(tǒng)計的Labeyrie法[8]，Knox-Thompson(簡稱K-T)法[9]和三階統(tǒng)計的斑點掩膜法[10]；空域重建法包括簡單位移疊加法、迭代位移疊加法以及基于選幀的相關位移疊加法等。本文采用K-T法實現(xiàn)對太陽圖像的高分辨率重建。

2.1　Labeyrie法

在滿足等暈性假設下，每幀短曝光圖像i(x)可以看作是目標o(x)和光學點擴展函數(shù)h(x)的卷積：

i(x)=o(x)*h(x)，

(1)

它的傅里葉變換可以表示為

I(f→)=O(f→)·H(f→)，

(2)

其中，I(f→)，O(f→)，H(f→)分別為i(x)，o(x)，h(x)的傅里葉變換；f為二維空間頻率變量。

Labeyrie法也稱為斑點干涉術，是通過統(tǒng)計斑點圖的平均能譜復原目標傅里葉振幅。按照Labeyrie的統(tǒng)計方法，多幀短曝光圖像頻域幅值平方的系綜平均為

(3)

2.2　K-T法

為了重建目標相位，Knox和Thompson于1974年對Labeyrie法改進后提出了一種基于計算頻譜的互相關獲取目標相位的方法，稱為Knox-Thompson(簡稱K-T)法，又稱互譜法。它在統(tǒng)計二階距時加上一個頻率的平移，則函數(shù)的互譜可以表示為

C(f→,Δf→)=I(f→)I*(f→+Δf→)，

(4)

其中，*表示共軛；Δf→表示空間頻率平移量，當Δf→=0時，互譜等同于能譜。對于Δf→的大小，應

斑點圖的互譜統(tǒng)計結果為

〈I(f→)I*(f→+Δf→)〉=O(f→)O*(f→+Δf→)·H(f→)H*(f→+Δf→)，

(5)

其中，H(f→)H*(f→+ Δf→)為K-T法的傳遞函數(shù)，它在接近衍射極限截止頻率時是一個顯著值的實數(shù)，所以互相關譜的相位代表目標各相鄰頻率處的相位差。通過從零頻處的迭代計算其他頻率處的相位值。

3　基于MPI-CUDA的K-T圖像重建算法實現(xiàn)

整個系統(tǒng)在單機環(huán)境下實現(xiàn)，計算流程如圖1，主要由3部分組成：MPI主進程，若干MPI子進程和圖形處理器程序。其中主進程負責分發(fā)和匯總數(shù)據(jù)，根據(jù)振幅和相位重建子圖并拼接成完整的視場圖。每個子進程負責計算若干子塊的振幅和相位，其中計算過程主要在圖形處理器上完成。

圖1MPI-CUDA實現(xiàn)高分辨重建流程圖
Fig.1Flowchartof MPI-CUDA implementation high resolution reconstruction

3.1　圖像預處理和分塊

在開始計算之前，首先要對圖像數(shù)據(jù)進行預處理，如平場和暗場修正。在成像系統(tǒng)中，線性空不變很難在一個完整視場內成立，但是在一個等暈區(qū)大小區(qū)域內線性空不變近似成立，可用相同的點擴展函數(shù)。因此可以根據(jù)等暈區(qū)把一組短曝光圖像切分成若干個子塊，每個子塊可以獨立地計算重建振幅和重建相位。最后拼接時只使用子塊的中間部分，因此切分時子塊之間有重疊，本文采用重疊50%。

3.2　數(shù)據(jù)分發(fā)

正是因為子塊之間的計算是不相關的，可以將每個子塊的重建看做獨立子任務。采用MPI并行計算框架，實現(xiàn)子塊重建的并行化，將所有子塊均衡地分到每個子進程中。為了減少MPI進程間的通信，尤其是主進程向各個子進程分發(fā)數(shù)據(jù)的時間消耗，將圖像原始數(shù)據(jù)保存在系統(tǒng)共享內存中。主進程只需向子進程傳遞一些控制信息，如子塊數(shù)據(jù)在圖像數(shù)據(jù)的索引位置等，每個子進程根據(jù)索引位置直接從共享內存中讀取圖像數(shù)據(jù)，從而節(jié)省大量通信時間。

雖然每個進程只負責子塊的計算，但是子塊的計算依然是數(shù)據(jù)密集的復雜計算。為了進一步對計算進行優(yōu)化，利用中央處理器-圖形處理器異構計算框架，在圖形處理器中實現(xiàn)數(shù)據(jù)的并行處理。在開始計算前，需要將數(shù)據(jù)拷貝到圖形處理器的顯存中，由于主機和圖形處理器之間通過PCI總線進行數(shù)據(jù)傳輸，每個子進程需循環(huán)處理多個子塊，頻繁的進行主機和設備之間數(shù)據(jù)傳輸嚴重影響性能。為了減小對計算性能的影響，在計算開始前把所有子塊的數(shù)據(jù)拷貝到顯存。

3.3　振幅重建

振幅的重建通常在頻域使用Labeyrie法退卷積傳遞函數(shù)完成，每個子塊的計算過程如下：

(1)計算斑點圖子塊的快速傅里葉變換；

(3)通過匹配求得相干長度r0；

首先需要將子塊圖像從空域轉換到頻域，使用CUDA的標準庫CUFFT對子塊進行傅里葉變換。步

3.4　相位重建

重建完整的圖像只有振幅是不夠的，還需要目標圖像的相位。采用K-T法計算得到圖像的互相關譜的相位即目標圖像的各相鄰頻域處的相位差，通過從零頻處的迭代得到其他頻率處的相位值。每個子塊的計算過程如下：

(1)計算斑點圖子塊的快速傅里葉變換；

(2)計算K個不同頻移的交叉譜；

(3)迭代計算K個不同頻移下的子塊相位；

(4)對K個相位結果求平均得到最終相位。

為了提高計算的精確度，對圖像進行K個不同的頻移并計算它們的交叉譜，由此帶來了計算量的大幅增加，采用串行的方式對每個頻率點進行計算則需要K×N個循環(huán)，可見計算量是巨大的。通過分析，不同頻移的交叉譜計算是獨立的，將線程數(shù)設置為與頻率點總數(shù)相等，每個線程計算一個頻率點的K個不同頻移。由互相關譜只能得到相鄰頻率間的相位差，要獲得目標相位還需要從零頻處迭代求和才能得到。頻率(m,n)處的相位值可以由不同路徑求得，如分別由頻率(m-1,n)和頻率(m,n-1)的相位與相對應的相位差求和得到，為了提高精度，通常對不同路徑的結果求平均。這種數(shù)據(jù)強依賴型的計算，每個點的相位都依賴前面的點計算得到，尚無較好的方法使用圖形處理器并行計算。

3.5　圖像拼接

根據(jù)得到各個子塊的目標傅里葉振幅和目標相位，通過逆傅里葉變換得到重建子圖像。重建完子圖像后，需要將所有子圖像拼接成全視場的圖像。因為使用頻域方法重建相位時，用的是疊加圖的初始相位，所以拼接時一般不會有偏移。拼接完成后，需要進行一些后續(xù)處理，如圖像的對齊等。

4　實驗分析

對上述CUDA-MPI的混合模型實現(xiàn)太陽圖像高分辨重建過程的性能進行測試。實驗采用的計算機配置為 IntelXeon-2620v3處理器，主頻為2.4 GHz；圖形處理器為兩張英偉達公司的TeslaK80顯卡，每張顯卡有兩個核芯，4 992個CUDA核，12 G顯存；操作系統(tǒng)是Ubuntu14.04.5LTS。測試數(shù)據(jù)為1 m太陽望遠鏡觀測的100幀短曝光太陽圖像，參數(shù)見表1。

表1　圖像數(shù)據(jù)參數(shù)Table 1　The parameters of the image

實驗中使用不同的進程數(shù)，測試MPI-CUDA和純MPI兩種并行方式實現(xiàn)K-T算法的時間，由于使用的顯卡共有4個計算核芯，為了使每個核芯的計算任務均衡，子進程數(shù)設置為4的倍數(shù)，結果見表2。從結果對比可以看出，只使用一個子進程即串行計算時，MPI-CUDA使得計算速度獲得明顯提升，整個流程的計算時間只有后者的1/3，這表明使用圖形處理器優(yōu)化計算后加速效果明顯。隨著子進程數(shù)目的增加，并行計算的優(yōu)勢展現(xiàn)出來，兩者的計算時間都大幅減少，但是MPI-CUDA的計算時間始終快于純MPI。不過受限于圖形處理器和中央處理器的計

表2　不同子進程數(shù)目下所用的時間Table 2　The time spent under the difference number of subprocesses

算能力，兩者的計算時間在達到一個最小值后，繼續(xù)增加進程數(shù)時間反而增加，其中MPI-CUDA在使用8個子進程時計算時間最少只有12 s，而MPI最少時卻需要20 s，因此MPI-CUDA混合并行的方式具有明顯優(yōu)勢。

此外在使用8個子進程的前提下，分別用中央處理器和圖形處理器處理子塊，并對各個環(huán)節(jié)所用的時間進行統(tǒng)計，結果如表3。

表3　中央處理器和圖形處理器計算時間對比Table 3　Computation time comparison between CPU and GPU

從時間對比上可以看出，算法移植到圖形處理器上后每一部分都獲得了速度提升，但是各部分的加速比卻有所差異，這和各部分在程序中的具體實現(xiàn)有關，并行化程度越高加速比越大，其中計算平均值獲得了172.56倍的加速。從時間占用上看，由于要計算k個不同頻移的交叉譜，交叉譜的計算時間占比最高，而使用圖形處理器加速后獲得了5.76的加速比，在單個子塊計算過程中的時間占用也從75.3%下降到了31.1%，對整個程序的加速起了關鍵作用。由于遞推求相位等數(shù)據(jù)依賴性強，計算無法在圖形處理器上并行實現(xiàn)，并且中央處理器還要完成圖像拼接等計算流程和邏輯控制流程，此外中央處理器和圖形處理器之間的數(shù)據(jù)傳輸也存在時間消耗，最終使用MPI-CUDA的方式較只使用MPI并行時間縮短了一半。由于圖形處理器與中央處理器結構的差異，大量的浮點數(shù)計算使得最后結果存在微小的精度誤差，但是MPI-CUDA混合并行模式得到的重建圖像達到了預期效果，圖像的紋理細節(jié)得以展現(xiàn)。

5　總　結

本文在單機上使用MPI-CUDA混合并行的模式實現(xiàn)了基于K-T算法的太陽高分辨圖像重建，并且取得了預期效果。實驗使用不同數(shù)目的進程進行測試，結果表明，這種混合并行的模式快于純MPI的并行，尤其是子塊中交叉譜的計算移植到圖形處理器后獲得了2.4倍的加速，使整個流程的時間大大減少，圖形處理器的密集運算能力得到了充分體現(xiàn)。天文觀測產生的數(shù)據(jù)量大，并且數(shù)據(jù)處理的算法復雜，中央處理器的計算能力遠遠不能滿足海量數(shù)據(jù)處理的需要，而中央處理器-圖形處理器這種高性能混合架構為天文研究的大規(guī)模計算任務提供了借鑒和參考。在未來的工作中，需要進一步對算法進行分析和優(yōu)化，如在不影響重建效果的情況下，減小子塊的重疊面積，從而減小子塊的數(shù)量。加深對圖形處理器計算架構的理解，將更多的計算在圖形處理器上完成，并將MPI-CUDA混合并行的方法推廣到采用斑點掩膜法的太陽高分辨圖像重建中。