鄧 濤，陳 東，代紅兵，王新華,3

(1. 云南大學(xué)信息學(xué)院，云南 昆明 650504；2. 中國科學(xué)院云南天文臺，云南 昆明 650216；3. 中國科學(xué)院大學(xué)，北京 100049)

太陽是影響人類活動最大的恒星，尤其是太陽活動對地球環(huán)境、氣候和天氣等的影響[1]。由于受到大氣湍流的影響，太陽光線在穿過地球大氣層時產(chǎn)生波前誤差，光路發(fā)生偏轉(zhuǎn)，觀測的太陽圖像出現(xiàn)不同程度的偏移、抖動、模糊等，導(dǎo)致地基望遠(yuǎn)鏡的觀測圖像質(zhì)量下降[2]。為了消除大氣湍流對望遠(yuǎn)鏡成像結(jié)果的影響，通常采用空間望遠(yuǎn)鏡、自適應(yīng)光學(xué)以及圖像重建技術(shù)等方式獲取太陽高分辨圖像。

目前圖像重建技術(shù)常用的算法主要有相位差法、多幀盲反卷積法、斑點干涉術(shù)、K-T算法、斑點掩模法、簡單位移疊加法、迭代位移疊加法以及選幀位移疊加法等。這些算法都是通過大量的短曝光圖像重建太陽高分辨圖像，往往因為數(shù)據(jù)量大、算法復(fù)雜導(dǎo)致無法達到實時重建的需求。

近年，國內(nèi)外在太陽高分辨圖像重建算法并行化研究方面做了很多工作，并獲得一定的加速比。文[3]基于雙節(jié)點集群，提出一種幀選擇和斑點掩模法的并行計算方法，采用多線程方法在57 s內(nèi)重建一幅256 × 256 pixel圖像。文[4]利用改寫于KISIP(Kiepenheuer Institut für Sonnenphysik)的算法，采用圖形處理器一次處理單組分塊數(shù)據(jù)，實現(xiàn)4.2 s內(nèi)重建225個128 × 128 pixel子塊的相位。文[5]從1 m新真空太陽望遠(yuǎn)鏡(New Vacuum Solar Telescope, NVST)的TiO通道選取一組子塊圖像(100 × 256 × 256 pixel)，采用OpenMP方法實現(xiàn)了一組子塊圖像的并行化，重建單幀256 × 256 pixel的子塊圖像，運行時間減少至2.7 s，獲得2.5倍加速。文[6]采用圖形處理器的統(tǒng)一計算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)對光球TiO通道中的一組子塊(100 × 256 × 256 pixel)實現(xiàn)并行化，基于統(tǒng)一計算設(shè)備架構(gòu)方法重建單組子圖的運行時間減少到約0.7 s。文[7]基于統(tǒng)一計算設(shè)備架構(gòu)在斑點掩模算法中實現(xiàn)單個子塊圖形處理器內(nèi)的并行化，采用并行方法比純中央處理器運行的串行算法加速比達到7。

綜上所述，斑點掩模法重建太陽高分辨圖像時的圖像并行化研究已經(jīng)取得一定的成果，但大多數(shù)局限于單純的中央處理器或圖形處理器加速，而且圖形處理器一次只能處理一組分塊數(shù)據(jù)，沒有完全發(fā)揮中央處理器/圖形處理器的并行化能力。如何將多組分塊數(shù)據(jù)分配到圖形處理器同時并行處理，進一步提高中央處理器和圖形處理器的并行計算能力和資源效率，本文提出了基于多進程并行加速的太陽高分辨圖像重建方法。

1 太陽高分辨圖像重建方法

1 m新真空太陽望遠(yuǎn)鏡[8]坐落在撫仙湖畔，主要的觀測波段有G-band(430.5 nm)、Hα(656.28 nm)和TiO(705.8 nm)。1 m新真空太陽望遠(yuǎn)鏡太陽高分辨圖像重建使用了兩個層次：Level1位移疊加法重建色球圖像，Level1+斑點掩模法重建光球或色球圖像。Level1+計算復(fù)雜度比Level1高，但Level1+在視寧度好的時候重建的圖像質(zhì)量更好，信噪比更高。

1 m新真空太陽望遠(yuǎn)鏡Level1+太陽高分辨圖像重建流程有圖像預(yù)處理、圖像初對齊、視寧度估計、圖像分塊處理和子塊拼接。其中，圖像分塊處理主要采用斑點干涉術(shù)和斑點掩模法重建太陽高分辨圖像的模和相位。圖像分塊處理包含若干環(huán)節(jié)，其中相位遞推環(huán)節(jié)數(shù)據(jù)量大，計算復(fù)雜，并且需要考慮多個遞推路徑的整合，在Level1+的重建流程中最為耗時。因此，子塊處理的好壞嚴(yán)重影響重建的效果和時效。

1.1 振幅重建

在滿足等暈區(qū)條件下，短曝光像是目標(biāo)和系統(tǒng)的點擴散函數(shù)的卷積

i(x,y)=o(x,y)?p(x,y)，

(1)

其中，?為卷積符號，時域中的卷積對應(yīng)于頻域中的乘積。(1)式在頻域中滿足

I(u,v)=O(u,v)P(u,v)，

(2)

其中，I(u,v)，O(u,v)和P(u,v)分別為時域中對應(yīng)項的傅里葉變換。

表5列出了2017年廣東省、江蘇省、北京市與上海市自然科學(xué)基金資助SCI論文合作發(fā)文量排名前10的國家和地區(qū)。廣東省與江蘇省、北京市、上海市自然科學(xué)基金資助SCI論文排名前10的合作國家和地區(qū)分布有相似之處，如與美國、澳大利亞、英國、加拿大、新加坡、日本、法國等經(jīng)濟、科技發(fā)達國家的合作較多。同時，4個省市自然科學(xué)基金資助SCI論文數(shù)排名前10的合作國家和地區(qū)也有各自的特點，如：除美國、澳大利亞、英國、加拿大、新加坡、日本、法國外，廣東省與中國臺灣地區(qū)和愛爾蘭合作較多；江蘇省與韓國和沙特阿拉伯王國合作較多；北京市與中國臺灣地區(qū)和新西蘭合作較多；上海市與韓國和比利時合作較多。

功率譜統(tǒng)計為

(3)

1.2 相位重建

目標(biāo)斑點圖的三重相關(guān)為

(4)

其中，ik(x),im(x)和in(x)為3個強度分布。目標(biāo)斑點圖的三重自相關(guān)的傅里葉變換

I(3)(u,v)=I(u)I(v)I(-u,-v)，

(5)

其中，I(u)是i(x)的傅里葉變換；u和v是二維空間頻率。

目標(biāo)斑點圖的平均重譜

〈I(3)(u,v)〉=O(3)(u,v)〈P(3)(u,v)〉，

(6)

其中，〈〉為系綜平均；〈P(3)(u,v)〉為平均斑點掩模法傳遞函數(shù)。在得到平均重譜后，由低頻到高頻的相位元逐步遞推，恢復(fù)目標(biāo)斑點圖的全部相位。

2 多進程并行加速方法

2.1 中央處理器/圖形處理器混合計算方法

為了滿足整個視場線性空間平移不變性，我們需要把預(yù)處理和初對齊后的圖像分割成一個個子塊，然后相同位置的子塊合并成子塊組。1 m新真空太陽望遠(yuǎn)鏡Level1+的太陽高分辨圖像重建中現(xiàn)有的并行加速方法大多局限于中央處理器/圖形處理器混合計算方法，即圖形處理器一次只能處理單個子塊組，不同子塊組之間還是串行計算?；谥醒胩幚砥?圖形處理器混合計算方法的圖像分塊處理流程如圖1。

雖然以上方法比純中央處理器圖像分塊串行處理的時間少，但這種方法存在兩個問題：(1)子塊組間的處理依次進行(串行)，在圖形處理器處理時，中央處理器處在空閑狀態(tài)，利用率不高；(2)圖形處理器一次只處理一個子塊組，利用率不高。為此，本文提出了基于多進程圖像分塊處理并行加速的方法。

2.2 多進程并行加速方法

針對現(xiàn)有方法存在的問題，本文引入多進程，提高中央處理器的利用率。多進程可以使圖形處理器處理更多的子塊組，同時提高中央處理器和圖形處理器的并行化處理能力。待處理的多幀圖像分割成很多子塊后，子塊合并為若干子塊組，把所有的子塊組加入任務(wù)列表，每個進程順序選擇任務(wù)列表中的一個子塊組，多個進程同時操作圖形處理器并行處理多個子塊組。基于多進程并行加速方法的圖像分塊處理流程如圖2。

圖2 基于多進程并行加速方法的圖像分塊處理流程Fig.2 Image block processing flow based on multi-process parallel acceleration method

圖1 基于中央處理器/圖形處理器混合計算方法的圖像分塊處理流程

多進程并行加速方法的圖像分塊處理流程為：

(1)創(chuàng)建任務(wù)列表。創(chuàng)建并初始化列表，將圖像分塊后的子塊組加入任務(wù)列表。

(2)創(chuàng)建進程池。創(chuàng)建并初始化一個進程池，并在進程池中添加適當(dāng)?shù)淖舆M程數(shù)量。

(3)分配進程任務(wù)。由主進程依次把任務(wù)列表中的任務(wù)分配給處于空閑狀態(tài)的子進程。

(4)傳遞參數(shù)。主要傳遞子塊組處理過程需要的數(shù)據(jù)和參數(shù)。

(5)處理子塊組。進程池開啟多少子進程，就有多少個子塊組同時并行計算。子塊對齊、模的重建、平均重譜的計算以及初始相位的計算等都在圖形處理器并行完成。圖形處理器完成初始相位計算后，把初始相位數(shù)據(jù)傳遞給中央處理器的相應(yīng)子進程開始相位遞推，完成相位遞推后，子進程將相位數(shù)據(jù)從中央處理器傳遞到圖形處理器，最后進行模和相位的合成。

(6)再次分配任務(wù)。處理完子塊組后處于空閑狀態(tài)的子進程都返回進程池，進程池將空閑子進程信息反饋給主進程，若子任務(wù)列表仍有待執(zhí)行的子塊組，主進程將待執(zhí)行的子塊組分配給進程池中空閑的子進程。

(7)關(guān)閉進程池。主進程反復(fù)檢測子塊組任務(wù)列表，當(dāng)任務(wù)列表中沒有待處理子塊組時，等待所有子進程執(zhí)行完畢，進程池將所有空閑子進程信息反饋給主進程，釋放所有子進程并關(guān)閉進程池。

由上述過程可知，各個子塊組相互獨立，互不影響，進程池中每個子進程處理相應(yīng)的子塊組，同時并行運行。子塊對齊、模的重建、初始相位計算以及模和相位的合成都是由多個子進程在圖形處理器并行完成的，但每個子塊組的相位遞推是由相應(yīng)的子進程在中央處理器并行計算的。值得注意的是，在創(chuàng)建和初始化進程池時，進程并不是越多越好，創(chuàng)建新進程會耗費系統(tǒng)資源，所以進程池中只能添加合適數(shù)量的子進程，數(shù)量取決于中央處理器和圖形處理器的資源能力。此外，為避免不必要的系統(tǒng)開銷，程序并不立即撤回完成任務(wù)的子進程，而是在任務(wù)列表還有任務(wù)時選擇進程復(fù)用，再次分配任務(wù)。

以Python為例的實現(xiàn)步驟如下：(1)多進程管理Multiprocessing模塊提供Pool類，通過import命令導(dǎo)入Multiprocessing模塊，使用multiprocessing.Pool函數(shù)導(dǎo)入Pool類。(2)創(chuàng)建并初始化任務(wù)列表blocklist=[]，將子塊組加入任務(wù)列表。(3)設(shè)置指定的工作子進程數(shù)目(num)，使用pool=Pool(num)，按照指定數(shù)量子進程創(chuàng)建和初始化進程池，供用戶調(diào)用。(4)使用pool.apply_async()，pool.map()，pool.apply()和pool.map_async()等方法，調(diào)用子塊組處理函數(shù)并傳遞子塊組任務(wù)列表參數(shù)，然后提交給進程池。當(dāng)新子塊組任務(wù)請求提交到進程池時，如果Pool不滿，則創(chuàng)建一個新進程執(zhí)行該任務(wù)請求，如果Pool已滿，告知進程池請求等待。(5)通過Cupy模塊中cupy.asarray()方法和cp.asnumpy()方法進行中央處理器和圖形處理器數(shù)據(jù)傳輸，每個子塊組計算完成后，由列表收集結(jié)果。(6)當(dāng)任務(wù)列表沒有待處理子塊組時，使用pool.close()方法，進程池不再接受新的任務(wù)；當(dāng)所有子塊組計算完成后，工作進程退出。(7)最后使用pool.join()方法，等待進程池中所有的子進程結(jié)束，返回主進程。關(guān)鍵過程的代碼實現(xiàn)如表1。

表1 關(guān)鍵過程的代碼實現(xiàn)Table 1 Code implementation of key processes

3 結(jié)果和分析

3.1 實驗環(huán)境

硬件：Intel(R)Core(TM)i5-8400 CPU@2.80 GHz(2 801 MHz)的處理器(6核)，8 GB隨機存取存儲器(Random Access Memory, RAM)，NVIDIA GeForce GTX 1050 Ti，4 095 MiB。軟件：Microsoft Windows 10，Spyder 4.0.13，CUDA 10.2，Python 3.7.6。

3.2 并行加速結(jié)果

以圖像分塊處理為例，我們對第2節(jié)兩種方法進行測時比較，實驗結(jié)果如表2。

表2 兩種方法處理時間Table 2 Processing time of two methods

為了驗證代碼的正確性，本文選取2020年6月6日1 m新真空太陽望遠(yuǎn)鏡Hα波段的觀測數(shù)據(jù)10組，每組100幀(1 028 × 1 024 pixel)，每幀采用重疊方式分割成96 × 96 pixel，將每幀相同位置對應(yīng)的子塊合并成一組，劃分為625個子塊組。我們還選取2020年6月8日1 m新真空太陽望遠(yuǎn)鏡TiO波段的觀測數(shù)據(jù)10組，由于受顯存大小的限制，每組選取50幀(2 160 × 2 560 pixel)，每幀采用重疊方式分割成128 × 128 pixel，每幀相同位置對應(yīng)的子塊合并成一組，劃分為1 050個子塊組。

實驗結(jié)果表明，Hα波段和TiO波段數(shù)據(jù)在采用多進程并行加速方法時，圖像分塊處理的平均時間相對較少，加速比分別為4.69和4.71。由于在中央處理器/圖形處理器混合計算方法中，子塊組之間仍然是串行計算，圖形處理器一次只能處理單組分塊數(shù)據(jù)；而在多進程并行加速方法中，子塊組之間是并行計算，圖形處理器能同時處理多組分塊數(shù)據(jù)?；诙噙M程并行加速方法可提高中央處理器和圖形處理器的資源利用效率，顯著提高圖像分塊處理的速度。

在多進程并行加速方法中，不同的進程數(shù)有不同的圖像分塊處理時間，如圖3。圖3顯示了不同的進程數(shù)對圖像分塊處理并行加速的耗時情況，隨著進程池中子進程數(shù)量的增加，Hα波段的所有子塊處理時間降低到157.76 s，TiO波段的所有子塊處理時間降低到203.32 s，并行加速的效果明顯。但是在使用6個進程以后，受到中央處理器數(shù)量、圖形處理器顯存以及多進程調(diào)度開銷的影響，圖像分塊處理時間上下浮動并呈現(xiàn)上升的趨勢。在選擇合適的子進程后，使用多進程并行加速方法處理時，中央處理器和圖形處理器利用率都提高了，中央處理器利用率可以達到100%。

圖3 圖像分塊處理在不同進程數(shù)下的耗時情況

4 結(jié) 語

針對現(xiàn)有方法子塊組間低效的串行處理導(dǎo)致中央處理器和圖形處理器利用率不高的問題，本文提出了基于多進程并行加速太陽高分辨圖像重建的方法，利用多核和多進程技術(shù)，有效提高了中央處理器/圖形處理器的利用率和重建速度，研究可以為天文數(shù)據(jù)并行化處理提供借鑒參考。進一步提高中央處理器/圖形處理器的并行化程度，仍然有許多亟待突破的關(guān)鍵問題，其中，中央處理器承擔(dān)的相位遞推壓力最大，耗費圖像分塊處理過程80%的時間。下一步我們考慮基于相位遞推的特點進行并行化，同時對相關(guān)算法進行優(yōu)化改進，使中央處理器/圖形處理器計算負(fù)載達到均衡。此外，相關(guān)研究還需要在消息傳遞接口(Message Passing Interface, MPI)/圖形處理器異構(gòu)環(huán)境中進行驗證。