Hα全日面云污染實(shí)時(shí)識(shí)別和修復(fù)系統(tǒng)*

張能維，楊云飛，李冉陽(yáng)，季凱帆

(昆明理工大學(xué)云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，云南 昆明　650500)

?

Hα全日面云污染實(shí)時(shí)識(shí)別和修復(fù)系統(tǒng)*

張能維，楊云飛，李冉陽(yáng)，季凱帆

(昆明理工大學(xué)云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，云南 昆明650500)

摘要：高空云層導(dǎo)致所觀測(cè)的Hα全日面像上覆蓋有一層云污染，使得圖像上的太陽(yáng)活動(dòng)細(xì)節(jié)變得模糊不清。為了能夠?qū)崟r(shí)探測(cè)云污染，并及時(shí)顯示修復(fù)后的圖像，采用圖形處理器技術(shù)實(shí)現(xiàn)了一個(gè)Hα全日面云污染實(shí)時(shí)識(shí)別和修復(fù)系統(tǒng)。該系統(tǒng)主要在統(tǒng)一計(jì)算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)環(huán)境下利用圖形處理器并行實(shí)現(xiàn)：(1)二值化圖像橢圓長(zhǎng)短軸比值法識(shí)別重度云污染圖像；(2)臨邊昏暗曲線中心對(duì)稱法識(shí)別可修復(fù)云污染圖像；(3)頻域巴特沃斯低通濾波法去除云污染。通過(guò)對(duì)系統(tǒng)中各運(yùn)算在圖形處理器中花費(fèi)的時(shí)間進(jìn)行詳細(xì)測(cè)量，發(fā)現(xiàn)傅里葉正反變換和頻域?yàn)V波占用了圖形處理器總處理時(shí)間的52.9%，是系統(tǒng)中最耗時(shí)的。然而，相對(duì)于1min的觀測(cè)時(shí)間間隔，約0.7 s的總處理時(shí)間可以滿足實(shí)時(shí)顯示的需要。另外，通過(guò)對(duì)修復(fù)后的圖像做質(zhì)量評(píng)價(jià)，驗(yàn)證了所采用的修復(fù)算法可以有效地去除云污染，并對(duì)太陽(yáng)活動(dòng)細(xì)節(jié)影響較小。最后討論了系統(tǒng)存在的不足和需要進(jìn)行的改進(jìn)。

關(guān)鍵詞：云污染；并行計(jì)算；識(shí)別和修復(fù)；巴特沃斯低通濾波器

Hα全日面觀測(cè)為監(jiān)測(cè)太陽(yáng)色球?qū)踊顒?dòng)現(xiàn)象提供大量的數(shù)據(jù)。全球擁有眾多的Hα全日面觀測(cè)站，如全球日震網(wǎng)(Global Oscillation Network Group, GONG)擁有全球6個(gè)觀測(cè)點(diǎn)[1]，國(guó)內(nèi)的北京懷柔觀測(cè)基地等等。雖然各站點(diǎn)在選址時(shí)進(jìn)行了嚴(yán)格的篩選，但是實(shí)際觀測(cè)中會(huì)出現(xiàn)大量的多云天氣[2]。在這樣的天氣觀測(cè)，經(jīng)常出現(xiàn)太陽(yáng)局部被云覆蓋的情況，導(dǎo)致圖像細(xì)節(jié)模糊。圖1顯示了Cerro Tololo美洲際天文臺(tái)在2015年4月20日觀測(cè)的3張圖像，可以看出后兩張圖像存在不同程度的云污染現(xiàn)象。

云覆蓋是導(dǎo)致全日面像畸變的一種情況。為了檢測(cè)和校正各種畸變，前人提出了一些關(guān)于Hα全日面像的質(zhì)量評(píng)價(jià)和修復(fù)算法。如文[3]采用雙尺度濾波歸一化圖像強(qiáng)度；文[4]根據(jù)圖像臨邊昏暗曲線的峰度判斷圖像質(zhì)量；文[5]根據(jù)最小相關(guān)系數(shù)法判斷圖像是否有云污染；文[6-7]采用大尺度中值濾波修復(fù)畸變圖像和去除云污染，文[5]使用多尺度形態(tài)學(xué)濾波法修復(fù)云污染圖像等。這些方法主要為了在后續(xù)的太陽(yáng)物理研究中，能在觀測(cè)圖像上更好地切割、識(shí)別和測(cè)量各種太陽(yáng)活動(dòng)特征結(jié)構(gòu)，如耀斑、暗條等等。因此，并不涉及執(zhí)行效率和運(yùn)行速度的問(wèn)題。

如果能將這些方法用于實(shí)時(shí)探測(cè)是否有云污染，并實(shí)時(shí)顯示修復(fù)后的圖像，將非常有利于觀測(cè)的進(jìn)行和觀測(cè)質(zhì)量的判斷。但這必然要求相應(yīng)的優(yōu)化算法和計(jì)算機(jī)的快速處理能力。為此基于文[5]的處理方法，采用圖形處理器技術(shù)實(shí)現(xiàn)了一個(gè)實(shí)時(shí)全日面像云污染的識(shí)別和修復(fù)系統(tǒng)。在硬件環(huán)境下，處理一張?jiān)莆廴緢D像消耗的總時(shí)間約為0.7 s。根據(jù)圖形處理器技術(shù)的特點(diǎn)，采用頻域巴特沃斯濾波器，并使用圖像相似度算法對(duì)無(wú)云圖像和云污染修復(fù)后的圖像進(jìn)行局部相似度檢驗(yàn)。結(jié)果顯示，在滿足處理速度的同時(shí)，這一濾波算法可以在很大程度上保留太陽(yáng)的活動(dòng)細(xì)節(jié)，滿足觀測(cè)實(shí)時(shí)顯示的需要。

圖1(a)、(b)、(c) 3幅圖像為Cerro Tololo美洲際天文臺(tái)分別于當(dāng)?shù)貢r(shí)間2015年4月20日 (16∶07、16∶08、20∶17) 拍攝

1算法原理分析

對(duì)云污染圖像的處理包括兩步：云污染的檢測(cè)和云污染的修復(fù)。其中檢測(cè)包括：(1)判斷全日面像是否為一個(gè)圓，否則認(rèn)為是重度云污染的圖像而放棄修復(fù)；(2)判斷全日面像的臨邊昏暗曲線是否圓對(duì)稱，如果不對(duì)稱，說(shuō)明有云污染。云污染的修復(fù)采用文[5]提出的模板法，即通過(guò)和一個(gè)標(biāo)準(zhǔn)全日面像模板進(jìn)行對(duì)比，并通過(guò)濾波得到云層透過(guò)率，從而修復(fù)被污染的圖像。

1.1重度云污染圖像的判斷

正常情況下的全日面像其太陽(yáng)圓面與天空背景的強(qiáng)度值相差很大，可以用一個(gè)高于天空背景的閾值對(duì)其進(jìn)行二值化。當(dāng)云污染過(guò)重時(shí)，日面像就會(huì)殘缺。圖2為圖1(c)二值化后的圖像。

因此，可先將二值化圖像假定為一個(gè)橢圓，通過(guò)計(jì)算其長(zhǎng)軸與短軸的比值E判斷是否殘缺。對(duì)良好觀測(cè)條件下的GONG圖像統(tǒng)計(jì)表明，其E值為1.0010 ± 0.0015。而圖2的E值為1.938，可判定為重度污染的圖像。在實(shí)際系統(tǒng)中，凡E值高于1.1的都?xì)w為重度污染的圖像。

1.2無(wú)云污染圖像的判斷

對(duì)于二值化后為一個(gè)圓的圖像，多數(shù)沒(méi)有云污染。根據(jù)文[5]的方法，可通過(guò)判斷臨邊昏暗曲線是否中心對(duì)稱判斷是否有云污染存在。主要步驟包括：(1)將太陽(yáng)圖像中心化并將其從直角坐標(biāo)系轉(zhuǎn)到極坐標(biāo)系；(2)獲取極坐標(biāo)圖像等分4部分中值臨邊昏暗曲線；(3)4條中值臨邊昏暗曲線兩兩求相關(guān)系數(shù)，比較得到最小相關(guān)系數(shù)。

圖2圖1(c)二值化處理后圖像

Fig.2Binary image processing result of Fig.1(c)

如果無(wú)云污染，這個(gè)系數(shù)接近1，反映了臨邊昏暗曲線的圓對(duì)稱性。當(dāng)有云污染時(shí)，這個(gè)系數(shù)會(huì)變小。在系統(tǒng)中采用的閾值為0.95。圖3顯示了圖1(b)的4條臨邊昏暗曲線，其最小相關(guān)系數(shù)為0.32，而相應(yīng)的圖1(a)為0.99。據(jù)此判定圖1中(a)不存在云污染，(b)存在可修復(fù)的云污染。

1.3云污染的修復(fù)

云污染全日面像的成像模型可簡(jiǎn)化為以下公式：

其中，c(x,y)為云層透射率；t(x,y)為沒(méi)有云污染的全日面像；T(x,y)為實(shí)際觀測(cè)的云污染全日面像。修復(fù)原理可簡(jiǎn)單表示成t(x,y)=T(x,y)/c(x,y)。修復(fù)的步驟為：(1)選擇一張無(wú)云圖像(如圖1(a))處理得到其標(biāo)準(zhǔn)臨邊昏暗輪廓(圖4(a))作為模板；(2)用有云污染的圖像(圖1(b))除標(biāo)準(zhǔn)臨邊昏暗輪廓模板，得到含有太陽(yáng)活動(dòng)細(xì)節(jié)的云層透射率(圖4(b))；(3)將圖4(b)轉(zhuǎn)換到頻域并利用巴特沃斯濾波器對(duì)其進(jìn)行低通濾波，扣除太陽(yáng)活動(dòng)細(xì)節(jié)的影響，得到云層透射率(圖4(c))；(4)將原始的云污染圖像除以云層透射率得到最終去除云污染的圖像(圖4(d))。

圖3圖1(b)對(duì)應(yīng)的四個(gè)象限的中值臨邊昏暗曲線

Fig.3The median limb-darkening profiles of four quadrants in azimuth of Fig.1(b)

圖4(a)由圖1(a)所得標(biāo)準(zhǔn)臨邊昏暗輪廓； (b)含太陽(yáng)活動(dòng)細(xì)節(jié)的云層透射率；(c)去除太陽(yáng)活動(dòng)細(xì)節(jié)的云層透射率； (d)去除云污染的圖像

Fig.4(a) The standard profile of Fig.1(a); (b) The clouds transmittance with solar features;(c) The clouds transmittance without solar features; (d) The solar image without cloud pollution

在第(3)步發(fā)現(xiàn)，圖4(b)中的太陽(yáng)活動(dòng)細(xì)節(jié)在頻域圖中主要分布于以中心為原點(diǎn)，半徑為D0的圓外區(qū)域。因此在巴特沃斯濾波器中，只需選擇合適的半徑D0(截止頻率)將其圓外的高頻部分濾除，即可在濾波后得到扣除太陽(yáng)活動(dòng)細(xì)節(jié)的云層透射率。經(jīng)過(guò)統(tǒng)計(jì)不同時(shí)段的云污染圖像，將系統(tǒng)中截止頻率取為圖像尺寸的0.003倍。

為充分發(fā)揮圖形處理器的并行處理優(yōu)勢(shì)，采用頻域?yàn)V波法。為此，采用圖像結(jié)構(gòu)相似度(Structural Similarity Index Measurement, SSIM)[8]算法證實(shí)修復(fù)的可靠性。選兩幅拍攝時(shí)間僅隔一分鐘的圖像，第1幅是干凈的圖像，第2幅為云污染的圖像，對(duì)第2幅圖像進(jìn)行修復(fù)得到第3幅圖像。圖5顯示了這3幅圖像的局部區(qū)域。計(jì)算得到無(wú)云圖像(圖5(a))與有云圖像(圖5(b))的圖像結(jié)構(gòu)相似度為0.299 7，無(wú)云圖像(圖5(a))與去除云圖像(圖5(c))的圖像結(jié)構(gòu)相似度為0.755 8?？梢钥闯鋈コ莆廴局髨D像的結(jié)構(gòu)相似性明顯高于去除云污染之前的圖像，說(shuō)明頻域?yàn)V波也可以相當(dāng)程度地去除云污染，并保留了太陽(yáng)的活動(dòng)細(xì)節(jié)。

2云污染圖像識(shí)別和恢復(fù)系統(tǒng)的圖形處理器實(shí)現(xiàn)

圖形處理器的設(shè)計(jì)初衷是為了將計(jì)算機(jī)中關(guān)于圖形圖像的處理和顯示從中央處理器中獨(dú)立出來(lái)。由于圖形處理器采用了單指令多數(shù)據(jù)流[9]的數(shù)據(jù)操作方式，因而并行計(jì)算能力很強(qiáng)。為了將圖形處理器的多線程并行能力運(yùn)用到通用計(jì)算上來(lái)，NVIDIA公司推出了統(tǒng)一計(jì)算設(shè)備架構(gòu)這一編程模型[10]。全日面像云檢測(cè)和修復(fù)算法具備數(shù)據(jù)計(jì)算密度高、可并行性強(qiáng)的條件，因而十分適合用圖形處理器處理。并且CUDA中有很多現(xiàn)成的函數(shù)庫(kù)，如CUBLAS庫(kù)中的大量針對(duì)矩陣的基本運(yùn)算函數(shù)以及CUFFT庫(kù)中的傅里葉正反變換函數(shù)，可以極大提高程序的運(yùn)行效率。例如，文[11]利用圖形處理器技術(shù)在國(guó)家天文臺(tái)懷柔觀測(cè)站實(shí)現(xiàn)一個(gè)實(shí)時(shí)深積分磁場(chǎng)觀測(cè)系統(tǒng)，文[12]利用圖形處理器技術(shù)實(shí)現(xiàn)了撫仙湖1m真空太陽(yáng)望遠(yuǎn)鏡level 1級(jí)圖像選幀。

圖5(a)、 (b)、 (c)分別為圖1(a)、 圖1(b)、 圖3(d)自中心截取的1 200 × 1 200像素局部圖

Fig.5(a), (b) and (c) are three 1200 × 1200 partial image, which are intercepted from the center of Fig.1(a), Fig.1(b) and Fig.3(d)

圖6展示了系統(tǒng)流程圖。該系統(tǒng)首先要顯示實(shí)時(shí)處理界面，然后循環(huán)檢測(cè)是否有更新的圖像。如果發(fā)現(xiàn)更新圖像先顯示原圖，再做云圖像檢測(cè)和修復(fù)處理，最后把處理后的圖像顯示在界面上與原圖像對(duì)比。如果沒(méi)有檢測(cè)到更新圖像，則界面一直顯示之前圖像的處理結(jié)果。

其中中央處理器所要完成的主要任務(wù)包括：(1)將原始FITS格式圖像數(shù)據(jù)從硬盤讀入中央處理器內(nèi)存并顯示在界面上；(2)將圖像數(shù)據(jù)從中央處理器內(nèi)存復(fù)制到圖形處理器顯存；(3)將處理后得到的數(shù)據(jù)從圖形處理器顯存復(fù)制回中央處理器內(nèi)存并顯示。

圖形處理器執(zhí)行系統(tǒng)的核心任務(wù)——云檢測(cè)和云去除。將圖形處理器中計(jì)算任務(wù)大致分為兩類：一類是通過(guò)直接調(diào)用CUDA庫(kù)函數(shù)實(shí)現(xiàn)。另一類是自行編寫(xiě)的利用線程并行實(shí)現(xiàn)的函數(shù)，其執(zhí)行的主要思想為：(1)根據(jù)數(shù)矩陣點(diǎn)的總數(shù)分配總線程數(shù)，通過(guò)線程號(hào)索引讓矩陣中每個(gè)點(diǎn)與線程一一對(duì)應(yīng)；(2)根據(jù)判斷條件將每個(gè)線程要完成的計(jì)算任務(wù)分類，同一判斷條件下所有線程的計(jì)算任務(wù)相同；(3)所有線程并行執(zhí)行各自分配的計(jì)算任務(wù)。

在檢測(cè)重度云污染圖像環(huán)節(jié)，首先并行實(shí)現(xiàn)圖像二值化處理。橢圓長(zhǎng)軸與短軸并行求取方法為：(1)并行求取橢圓重心位置；(2)并行得到橢圓上各點(diǎn)與其重心的協(xié)方差。最后將所得協(xié)方差元素代入公式即可得到橢圓長(zhǎng)短軸之比。據(jù)此如果判斷圖像云污染過(guò)重，則丟棄，否則繼續(xù)下面步驟。

接下來(lái)要檢測(cè)圖像是否存在可修復(fù)的云污染。首先并行實(shí)現(xiàn)將圖像上的全日面從日心(對(duì)應(yīng)上面所求重心)位置平移到圖像中心位置，即得到中心化圖像。中心化過(guò)后，再并行實(shí)現(xiàn)用雙線性插值法將圖像轉(zhuǎn)換至極坐標(biāo)系下。緊接著對(duì)極坐標(biāo)圖像角度的4個(gè)等間距區(qū)域分別進(jìn)行排序，排序方法是并行實(shí)現(xiàn)的雙調(diào)排序算法。從排完序的極坐標(biāo)圖像中直接取出所需的4條臨邊昏暗曲線，將兩條曲線用現(xiàn)有函數(shù)進(jìn)行點(diǎn)積、求和等運(yùn)算，再將結(jié)果代入公式即可得兩條曲線的相關(guān)系數(shù)。利用比較所得最小相關(guān)系數(shù)如判斷有云污染，則不需要修復(fù)，否則繼續(xù)下面去除云污染的步驟。

在使用模板法去除云污染時(shí)，要用到模板的標(biāo)準(zhǔn)臨邊昏暗輪廓。生成模板算法與檢測(cè)云算法很相似，圖7(a)展示了Hα全日面像生成模板系統(tǒng)的界面。先通過(guò)并行相除得到帶有太陽(yáng)活動(dòng)細(xì)節(jié)的云層透射率。緊接著利用CUFFT庫(kù)中函數(shù)進(jìn)行傅里葉變換將圖像轉(zhuǎn)換到頻域。再并行生成一個(gè)巴特沃斯低通濾波器，將其與頻域圖像并行相乘，即可得到濾除高頻部分的頻域圖像。然后將濾除高頻部分的頻域圖像利用傅里葉逆變換轉(zhuǎn)回時(shí)域得到消除太陽(yáng)活動(dòng)細(xì)節(jié)的云層透射率。最后將中心化云污染像與云層透射率并行相除得到最終去除云污染的圖像。

對(duì)于系統(tǒng)中一些重要參數(shù)，根據(jù)以往的數(shù)據(jù)實(shí)驗(yàn)統(tǒng)計(jì)設(shè)置默認(rèn)值。為了讓用戶可以在合理范圍內(nèi)進(jìn)行調(diào)整，同時(shí)將這些參數(shù)框放在界面上。如圖7(b)為Hα全日面像云污染實(shí)時(shí)識(shí)別和修復(fù)系統(tǒng)界面。

圖6 　Hα全日面云污染實(shí)時(shí)識(shí)別和修復(fù)系統(tǒng)

圖7　(a)生成模板系統(tǒng)界面；(b)云污染實(shí)時(shí)識(shí)別和修復(fù)系統(tǒng)界面

3系統(tǒng)性能測(cè)試

實(shí)驗(yàn)采用的計(jì)算機(jī)配置為：中央處理器為英特爾Pentium(奔騰)雙核E5200處理器，圖形處理器卡為Nvidia公司的GeForce GT 430。操作系統(tǒng)為Windows 7旗艦版64位。軟件開(kāi)發(fā)環(huán)境為Microsoft Visual Studio 2010下基于對(duì)話框的微軟基礎(chǔ)類庫(kù)(Microsoft Foundation Classes, MFC)和CUDA 6.5。

對(duì)系統(tǒng)中各主要運(yùn)算在圖形處理器中花費(fèi)的時(shí)間進(jìn)行測(cè)量，測(cè)試數(shù)據(jù)為圖1(b)云污染圖像(尺寸為2 048×2 048像素)，得到如表1的測(cè)試結(jié)果。

表1　系統(tǒng)中各運(yùn)算在圖形處理器中處理時(shí)間

經(jīng)測(cè)試， 除了圖形處理器的處理時(shí)間外，中央處理器的總運(yùn)行時(shí)間為342 ms。其中，中央處理器花費(fèi)的時(shí)間包括：從硬盤讀取原始圖像數(shù)據(jù)和模板圖像數(shù)據(jù)(154 ms)、在界面上顯示原始圖像和修復(fù)后的圖像(160 ms)、其余運(yùn)算(28 ms)。所以，該系統(tǒng)處理總時(shí)間為718 ms。

4討論

圖形處理器對(duì)高密度圖像數(shù)據(jù)的并行處理能力在系統(tǒng)中得到了充分體現(xiàn)。對(duì)于一幅2 048×2 048像素的云污染Hα全日面圖像，系統(tǒng)在圖形處理器中完成從檢測(cè)到修復(fù)等處理運(yùn)算的總時(shí)間為376.1ms。另外，系統(tǒng)在中央處理器中完成讀圖像、顯示圖像等處理總耗時(shí)為342 ms。針對(duì)常規(guī)1min采樣間隔的Hα全日面觀測(cè)，可以滿足觀測(cè)過(guò)程中對(duì)處理結(jié)果實(shí)時(shí)顯示的需要。為了發(fā)揮圖形處理器的并行優(yōu)勢(shì)，在修復(fù)過(guò)程中采用頻域巴特沃斯低通濾波的方法。通過(guò)使用圖像結(jié)構(gòu)相似度算法對(duì)修復(fù)后圖像做質(zhì)量評(píng)價(jià)，發(fā)現(xiàn)該方法與前人所用的中值濾波法和多尺度形態(tài)學(xué)濾波法修復(fù)效果相當(dāng)。系統(tǒng)除了可以處理GONG的數(shù)據(jù)，還使用懷柔太陽(yáng)觀測(cè)基地的Hα全日面像進(jìn)行測(cè)試，測(cè)試結(jié)果正常。此外，這個(gè)系統(tǒng)不但可以用于Hα全日面像的云污染處理，也可用在其他全日面像的處理中。

雖然統(tǒng)一計(jì)算設(shè)備架構(gòu)的Thrust函數(shù)庫(kù)提供了對(duì)于任意長(zhǎng)度數(shù)據(jù)的排序函數(shù),但需要先將數(shù)據(jù)放入向量容器，且無(wú)法實(shí)現(xiàn)數(shù)據(jù)集間的并行，導(dǎo)致排序時(shí)間過(guò)長(zhǎng)(單個(gè)數(shù)據(jù)集耗時(shí)0.5 ms，總耗時(shí)1.7 s)。而采用雙調(diào)排序不但數(shù)據(jù)集內(nèi)部的排序效率高于庫(kù)中的函數(shù)，而且可以實(shí)現(xiàn)數(shù)據(jù)集之間并行(受共享顯存的限制，單次可并行數(shù)據(jù)集有限，經(jīng)測(cè)試每次并行12個(gè)數(shù)據(jù)集總耗時(shí)最短，為26.8 ms)。雖然該算法對(duì)所排數(shù)據(jù)的長(zhǎng)度有約束(長(zhǎng)度必須為2的冪，因此將極坐標(biāo)角度設(shè)為512)，但并不影響后續(xù)的判斷。

在云污染的修復(fù)環(huán)節(jié)，生成巴特沃斯低通濾波器用時(shí)45.6 ms，相對(duì)于除傅里葉正反變換的其他運(yùn)算來(lái)說(shuō)，這是一個(gè)十分耗時(shí)的運(yùn)算。主要原因是每個(gè)線程所承擔(dān)的運(yùn)算量較大，而圖形處理器單個(gè)處理單元性能不及中央處理器內(nèi)核。曾嘗試將該濾波器做成一張圖像保存在硬盤中，讓系統(tǒng)每次從硬盤中讀取然后復(fù)制到圖形處理器中參與運(yùn)算。但測(cè)試發(fā)現(xiàn)，讀取用時(shí)77 ms，再?gòu)?fù)制到圖形處理器用時(shí)14.2 ms，總耗時(shí)約為在圖形處理器中生成的兩倍。

然而，系統(tǒng)也存在一些不足之處。首先是檢測(cè)可修復(fù)云污染算法，針對(duì)云污染對(duì)稱均勻分布的情況檢測(cè)失效(當(dāng)然這種情況極少出現(xiàn))。其次是頻域?yàn)V波的缺陷：(1)當(dāng)云污染比較嚴(yán)重時(shí)，如果通過(guò)增大濾波器的截止頻率更大限度地去除云污染，則同時(shí)會(huì)去除很多太陽(yáng)活動(dòng)細(xì)節(jié)：(2)圖像在太陽(yáng)邊緣處存在較大梯度差使得濾波時(shí)出現(xiàn)一些邊緣效應(yīng)。

要強(qiáng)調(diào)系統(tǒng)的主要目的是在觀測(cè)時(shí)實(shí)時(shí)顯示云污染修復(fù)后的太陽(yáng)圖像，因而著重關(guān)注處理速度。雖然對(duì)修正后的數(shù)據(jù)也進(jìn)行了保存，但如果需要進(jìn)一步的對(duì)這些云污染圖像進(jìn)行高精度處理，如對(duì)各種太陽(yáng)活動(dòng)現(xiàn)象進(jìn)行自動(dòng)識(shí)別、分割和測(cè)量，采用文[5]提出的多尺度形態(tài)學(xué)濾波更為合適。如果將多尺度形態(tài)學(xué)濾波并行化作為后續(xù)需要考慮的工作，同時(shí)也需要考慮對(duì)由于不精確的平場(chǎng)改正、快門效應(yīng)以及圖像飽和等原因所致畸變的Hα全日面圖像的修復(fù)，使得系統(tǒng)有更廣的適用范圍。進(jìn)一步驗(yàn)證系統(tǒng)涉及的云污染識(shí)別與修復(fù)原理是否同樣適用于太陽(yáng)局部像的云污染處理，如果適用再利用圖形處理器技術(shù)開(kāi)發(fā)相應(yīng)的處理軟件。

致謝：同時(shí)感謝由NOS主持的GONG項(xiàng)目和國(guó)家天文臺(tái)懷柔觀測(cè)基地所提供的Hα全日面圖像。

參考文獻(xiàn)：

[1]The GONG network of sites[EB/OL]. [2015-06-09]. http://gong.nso.edu/sites/.

[2]Hill F, Fischer G, Forgach S, et al. The Global Oscillation Network Group site survey, 2: results[J]. Solar Physics, 1994, 152(2): 351-379.

[3]Fuller N, Aboudarham J, Bentley R D. Filament recognition and image cleaning on Meudon Hα spectroheliograms [J]. Solar Physics, 2005, 227(1): 61-73.

[4]羅鷺鷥, 楊云飛, 季凱帆, 等. 基于臨邊昏暗輪廓的Hα全日面太陽(yáng)圖像質(zhì)量評(píng)價(jià)[J]. 天文學(xué)報(bào), 2013, 54(4): 392-400.

Luo Lusi, Yang Yunfei, Ji Kaifan, et al. The quality assessment of Hα full disk solar image based on limb darkening profile[J]. Acta Astronomica Sinica, 2013, 54(4): 392-400.

[5]Feng Song, Lin Jiaben, Yang Yunfei, et al. Automated detecting and removing cloud shadows in full-disk solar images[J]. New Astronomy, 2014, 32(10): 24-30.

[6]楊云飛, 馮松, 季凱帆. 畸變?nèi)彰嬗^測(cè)圖像的修復(fù)[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2014, 11(1): 54-59.

Yang Yunfei, Feng Song, Ji Kaifan. A new approach for correcting distortions in full-disk solar images[J]. Astronomical Research and Technology——Publications of National Astronomical Observatories of China, 2014, 11(1): 54-59.

[7]朱海波, 楊云飛, 鄧輝, 等. Hα全日面太陽(yáng)圖像云污染的去除[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2014, 11(2): 134-139.

Zhu Haibo, Yang Yunfei, Deng Hui, et al. Removal of cloud extinction for Hα full-disk solar image[J]. Astronomical Research and Technology——Publications of National Astronomical Observatories of China, 2014, 11(2): 134-139.

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

[9]Goyal N, Ormont J, Sankaralingam K, et al. January 2008 signature matching in network processing using SIMD/GPU architectures[R]. Cs Technical Reports, 2008.

[10]Kirk D. NVIDIA CUDA software and GPU parallel computing architecture[C]// Proceedings of the 6th International Symposium on Memory Management. 2007: 578.

[11]Shen Yangbin, Lin Jiaben, Ji Kaifan, et al. New real-time correlation solar observing system based on GPU for acquiring the deep-integration magnetogram[J]. New Astronomy, 2013, 25(2): 32-37.

[12]施正, 向永源, 鄧輝, 等. 一米真空太陽(yáng)望遠(yuǎn)鏡Level1級(jí)圖像選幀的GPU實(shí)現(xiàn)[J]. 科學(xué)通報(bào), 2015, 60(15): 1408-1413.Shi Zheng, Xiang Yongyuan, Deng Hui, et al. A method of Level 1 frames-selection based on GPU for new vacuum solar telescope[J].Chinese Science Bulletin, 2015, 60(15): 1408-1413.

A Real-Time Image Processing System for Detecting and Removing Cloud Shadows on Hα Full-Disk Solar

Zhang Nengwei, Yang Yunfei, Li Ranyang, Ji Kaifan

(Computer Technology Application Key Laboratory of Yunnan Province, Kunming University of Science and Technology, Kunming 650500, China, Email: flyingwithcloud@sina.com)

Abstract:Sky clouds lead to Hα full-disk solar observed images covered with shadows, which obscure the details of solar features. In this paper, a real-time image restoring system is introduced. For detecting and removing cloud shadows, and displaying the restoring image in real-time, Graphic Processing Unit (GPU) with Compute Unified Device Architecture (CUDA) parallel programming environment is employed to accelerate computing in the system. The parallel computing steps include: (1) identifying heavy cloud covered images by calculating the ratio of major and minor axes of a fitted ellipse; (2) recognizing recoverable cloud-covered images by calculating the symmetric of limb darkening curve; (3) removing cloud shadows by Butterworth low-pass filter in frequency domain. The total processing time is about 0.7s for one cloud-covered image; it meets the need of displaying in real-time. The processing time of each operation taken in GPU is measured. The results show that Fourier transform, inverse Fourier transform and Butterworth low-pass filtering are the most time-consuming part and together they take up 52.9% of the total time. This paper also evaluate the quality of restored image, and the conclusion is that the algorithm can effectively eliminate cloud shadows; meanwhile, it has little effect on the details of solar features. The existing problem and probable improvement are also discussed at the end of this paper.

Key words:Cloud shadows; Parallel computing; Detecting and removing; Butterworth low-pass filtering

基金項(xiàng)目：國(guó)家自然科學(xué)基金 (11303011, 11263004, 11463003, 11163004, U1231205) 資助.

收稿日期：2015-06-18；

修訂日期：2015-07-13

作者簡(jiǎn)介：張能維，男，碩士. 研究方向：天文圖像處理. Email: 281689743@qq.com 通訊作者：楊云飛，女，副教授. 研究方向：計(jì)算機(jī)應(yīng)用與天文圖像處理. Email: flyingwithcloud@sina.com

中圖分類號(hào)：P111; TP274

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào):1672-7673(2016)02-0242-08

CN 53-1189/PISSN 1672-7673