李洪波 趙明宇 劉煜1

(1中國科學院云南天文臺 昆明 650011) (2中國科學院大學 北京 100049)

基于相位一致性的冕環(huán)識別方法?

李洪波1,2?趙明宇1?劉煜1

(1中國科學院云南天文臺 昆明 650011) (2中國科學院大學 北京 100049)

嘗試將相位一致性方法增強的圖像用于冕環(huán)結構的識別.針對相位一致性方法增強圖像,提出一個新的冕環(huán)識別方法(簡稱PCB),該方法依據冕環(huán)形態(tài)變化比較平緩的形態(tài)學特性,將冕環(huán)識別方向限制在一個較小的范圍之內,并且首次同時將識別過程中推進方向的變化情況和識別點的相位一致性大小作為終止識別的判據.對多組數據的識別和分析則表明相位一致性方法增強的圖像的確適合用于冕環(huán)識別,而PCB方法所識別的冕環(huán)結構則同時具有很好的完整性和準確性,是一套現(xiàn)實可行的冕環(huán)自動識別方法.

太陽:日冕,技術:圖像處理,方法：觀測

1 引言

冕環(huán)作為日冕中主要的X射線發(fā)射源,被認為是由磁場束縛的熱等離子體形成的亮環(huán)狀結構[1].因此冕環(huán)的結構很好地反映了日冕的磁場結構,同時也是光球層和日冕之間的能量通道.冕環(huán)普遍存在于日冕層,可被大量的空間衛(wèi)星探測到,對冕環(huán)的研究有利于我們更好地理解日冕中等離子體的特性[2-3].關于冕環(huán)的基本性質和其演化規(guī)律卻一直未被很好地理解[4].首先各種冕環(huán)的分類標準還沒有很好地建立起來.初期Vaiana等[5]依據冕環(huán)的形態(tài)學特性為標準,將冕環(huán)分為亮點冕環(huán)、活動區(qū)冕環(huán)和大尺度冕環(huán).隨著多波段觀測數據的日益豐富,更詳細的溫度和密度研究得以展開.由于不同溫度冕環(huán)的溫度和密度性質明顯不同,例如暖冕上不同位置的溫度大致相同,而熱冕環(huán)則不同,另一種新的基于溫度的分類方法被提出,依據溫度的不同將冕環(huán)分為冷冕環(huán)[6]、暖冕環(huán)、熱冕環(huán)[7]3種.不同溫度冕環(huán)溫度和密度特性的不同也暗示著這些冕環(huán)可能產生于不同的物理機制[8].至于不同種類的冕環(huán)產生的具體物理機制以及它們是否可以相互轉化等基本問題,都未被很好地解決[9].另外,由于目前直接測量日冕磁場非常困難,對冕環(huán)幾何形狀的研究在探索日冕磁場結構和日冕中磁重聯(lián)方面也有著很重要的參考意義.大體上看,冕環(huán)呈現(xiàn)出對稱的環(huán)狀結構[10],比較特別的是冕環(huán)的橫截面沿著冕環(huán)有著很小的變化[11],這與我們通常使用的日冕中的無力場假設所預期的結果有些不同,無力場外推結果預測磁場結構的截面將會有幾倍的擴展.基于這個問題,有不同的冕環(huán)模型被提出,但是還沒有任何一個模型可以解釋所有的問題.值得注意的是,由于我們的觀測結果是一個投影效果,所以冕環(huán)的真實形狀的研究也有一定的難度,而“日地關系觀測臺”(Solar Terrestrial Relations Observatory,STEREO)[12]的發(fā)射為我們提供了多角度太陽觀測數據,這使得我們可以更好地重構出三維冕環(huán)結構[13].除此之外,冕環(huán)內部的等離子體流動[14–16],以及磁流體波動在冕環(huán)中的傳播等問題[17–19],都沒有被很好地理解.更加詳細的研究,以及更加深入的學習,毫無疑問都需要大量的樣本,這就需要我們對冕環(huán)進行大量準確的識別和提取.

自從冕環(huán)結構被發(fā)現(xiàn)以來,已經有多個團組對冕環(huán)結構的提取以及自動識別方法展開研究,并提出了幾種不同的冕環(huán)識別方法.主要的冕環(huán)自動識別方法基本上包括兩個主要的過程:首先要對原始圖像進行處理從而使冕環(huán)結構更加突出,或者對屬于冕環(huán)結構的點進行標記;第2步就要通過一些邊緣連接或者邊緣檢測程序對處理過的圖像中的冕環(huán)結構進行提取和識別.早在2006年Lee等[20]就首先提出了一種基于方向連通性(oriented connectivity)的自動冕環(huán)識別方法(簡稱OCM).這種方法首先將原始圖像進行銳化,然后用Strous提出的邊緣檢測方法1h ttp://www.lm sal.com/aschw and/stereo/2000easton/cdaw.h tm l對銳化圖像中的冕環(huán)結構進行標記,而接下來對標記的冕環(huán)結構進行識別連通時,作者建設性考慮了磁場的物理特性,并用一個磁場模型來指導接下來的識別過程.測試結果則表明OCM方法識別的冕環(huán)有著很好的完整性,識別的冕環(huán)也比較光滑.但是可能由于磁場模型本身和冕環(huán)實際情況之間的差別,以及第1步邊緣檢測的結果連續(xù)性較差而且噪點較多等因素的影響,盡管識別過程中兩次使用樣條擬合(sp line fitting),該程序所識別的冕環(huán)數目還是偏少.另外由于兩次人為擬合,也導致OCM方法的靈敏度變得很低,有時候甚至會將相鄰的冕環(huán)識別為一個冕環(huán).之后Aschwanden等[21]也提出一種冕環(huán)識別方法(簡稱ODM),該方法直接用冕環(huán)結構提升處理之后的圖像作為參考進行第2步的冕環(huán)結構識別,之后在沿著預測方向一定距離(假設5像素)附近一個較大的盒子(假設10 pixel×10 pixel)里面進行搜索,該預測方向是通過對每一個搜索盒子里最大值的位置進行線性回歸擬合(linear regression fitting)得到,將該窗口內最大亮度值以及窗口內最大亮度位置到該方向的距離兩個參數用作停止識別的判據.由于采用更直接的冕環(huán)增強圖像作為參考和新的識別方法,使得Aschwanden等的ODM冕環(huán)識別方法可以識別較多的冕環(huán)結構,但是測試結果也表明該方法也存在一定的不足,即抗噪聲能力較差,識別過程很容易受到噪聲影響.這使得其識別的冕環(huán)完整性較差,此外其結果也存在一些錯誤識別的小結構.

目前已經有多種針對邊緣結構提取的方法被開發(fā)出來,但是多數方法是針對一個更加廣泛的范圍,而專門針對冕環(huán)結構提取的方法則相對較少.因此在2008年Aschwanden等[21]將現(xiàn)有的5種對邊緣結構進行自動識別的方法進行了對比.除了前文介紹的兩個方法以外,他們還介紹了3種不同的結構識別方法.其中包括Lee等[20]開發(fā)的Dynam ic Aperture-based Loop Segmentation(簡稱DAM)方法,這種方法主要的特點是通過沿冕環(huán)的強度曲面是連續(xù)的高斯曲面這樣的假設來進行冕環(huán)的檢測和識別.該方法首先搜索圖像中符合高斯曲面的結構視為冕環(huán)結構,接著將結果中相鄰的結構依據擬合參數的相似情況進行連接,最后還要對不是冕環(huán)的結構進行剔除.這種DAM方法的結果則同前面介紹的OCM的結果類似,由于有連接的過程,所以識別的冕環(huán)會比較完整,但是其同樣存在靈敏度較低的問題,識別的冕環(huán)數目也比較少.文中另外介紹的一個冕環(huán)識別方法是Steger等2Steger L H.Techn ica l Rep ort FGBV-96-03,Forschungsgruppe B ildverstehen(FG BV),In form atik IX, Techn ische Un iversit¨at M¨unchen,G erm any,1996提出的一種邊緣探測方法(簡稱UDM),該方法已經被成功地用于遙感、醫(yī)學等多個領域.在這種方法中曲線結構的識別是基于該結構垂直方向的二階導數來實現(xiàn)的,這樣可以減小結構寬度以及不對稱性的影響,這種方法的結果則同Aschwanden的ODM方法類似,有著較好的靈敏度,可以識別出更多種類的冕環(huán),但是識別冕環(huán)的完整性則較差,其識別結果同樣存在一些奇怪的小尺度結構.除此以外Aschwanden等[21]還介紹了一種自動探測的方法(簡稱RAS)[22],這種方法是Inhester等人通過改進Parent等[23]的邊緣探測方法得到的,其主要是用圖像的泰勒級數(Taylor coefficient)來確定冕環(huán)點的位置、方向等參數.測試表明該方法不但可以識別不同種類的冕環(huán)結構,而且識別的冕環(huán)也比較完整,但是RAS的識別暗條的末端有時會存在一些錯誤彎曲,同時也會錯誤地識別一些小尺度的非冕環(huán)結構.Aschwanden的研究將這5種現(xiàn)存的不同冕環(huán)結構識別方法應用于“過渡區(qū)和日冕探測器”(Transition Region and Coronal Exp lorer,TRACE)[24]的171?A波段觀測圖像中,對提取的結果進行了詳細的比較.結果表明盡管5種算法都可以得到大致相同的結果,但同人工識別的結果相比,這些結果間還是存在著不可忽略的差異.因此如果想得到更好的結果,仍然需要更加詳細的研究[21].

值得注意的是,前面我們介紹的大量關于冕環(huán)自動識別的工作主要是優(yōu)化冕環(huán)自動識別工作的第2個步驟,針對冕環(huán)標記、邊緣連接和識別方法進行研究優(yōu)化,而對原圖像進行處理從而提高冕環(huán)結構對比度方面則沒有太多的方法改進.在這方面M cAteer等[25]測試了一種基于小波變換的冕環(huán)探測方法(簡稱W TMM),他們嘗試利用小波變換方法對冕環(huán)結構進行增強,并且針對TRACE在171?A波段的觀測圖像進行了研究.結果表明利用多尺度頻域信息增強后的冕環(huán)結構具有很好的方向連通性,可直接提取感興趣尺度信息等天然優(yōu)勢,但這種方法也存在著對較暗結構的提升效果較差、識別準確性較低等問題.類似于M cAteer等[25]的小波變換的方法,相位一致性同樣是基于圖像多尺度頻域信息的結構增強方法,不同的是該方法主要利用圖像的相位一致性反映圖像的結構信息、受到結構對比度影響小等特點來提取圖像的結構特征. Feng等[26]在2014年關于相位一致性方法在探測低對比度太陽結構方面的研究則很好地表明相位一致性方法在增強日冕邊緣冕環(huán)結構方面有著很大的優(yōu)勢.基于相位一致性方法的這些特性,我們發(fā)現(xiàn)相位一致性方法在冕環(huán)圖像處理方面的應用非常值得進一步研究,因此我們對相位一致性方法在冕環(huán)結構提取方面的應用進行了進一步的研究,并且提出了一個新的基于相位一致性的冕環(huán)識別方法.

2 相位一致性方法

研究表明圖像的結構特征更多地包含在圖像的相位譜之中,早在1981年Oppenheim和Lim[27]就做過相關的實驗.試驗中作者首先將兩幅圖像進行傅里葉變換分別計算出相位和振幅譜,然后將兩幅圖的功率譜進行交換并且反變換之后發(fā)現(xiàn)得到新的圖像與為其提供相位譜的圖像非常相似,而與為其提供功率譜的圖像有很少相似之處.

在結構識別的研究中,Morrone和Owens[28]將相位一致性定義為:

其中PC(x)代表x點的相位一致性,An(x)為x點第n個傅里葉分量的振幅,E(x)是x點上各個傅里葉分量的矢量和.它的幾何描述如圖1所示,一維信號F(x)在x處每一個傅里葉分量都有一個振幅An(x)和相位?n(x),這些傅里葉分量作為復矢量按照平行四邊形法則加起來,得到從原點到終點的復矢量振幅|E(x)|,相位一致性就是這個總振幅與各個傅里葉分量振幅之和的比值.顯然如果所有傅里葉分量相位角都相等,則PC(x)=1,而各個傅里葉分量相位完全不相干,則PC(x)趨近于0.

圖1 相位一致性的幾何描述[29]Fig.1 Geom etric descrip tion of the phase congruency[29]

值得注意的是,這里的相位一致性僅僅在一個較寬的頻域范圍內才有意義,如果某一個振幅分量遠遠大于其他振幅分量,則相位一致性結果總是趨近于1.因此Kovesi在相位一致性的計算公式中乘入了一個權重函數,從而得到了相位一致性新的表達式[29-30]:

這里??n(x)=cos(?n(x)??(x))?|sin[?n(x)??(x)]|,而T是估計的噪聲,ε為一個小量,為了防止分母為零的情況出現(xiàn),w(x)即新加入的權重因子,它的表達式為:

其中γ為增益因子,用來控制截止函數的陡峭程度,c則是一個截止參量,如果f(x)小于c則PC的值將被限制變小.N是濾波函數尺度的總個數,Amax(x)是各階傅里葉分量振幅的最大值.

上面主要描述的是一維信號的處理方式,而對于二維圖像的相位一致性,目前有兩種方法可以計算:一個就是方向濾波方法[29-30],另外一個則是單演濾波方法(Monogenic Filter)[31].這里我們采用的是Felsberg和Sommer提出的單演濾波方法[31].該方法在頻域中引入濾波函數:

式中i為虛數單位.為了減少計算量,這里引入濾波函數H,使H滿足

這里采用頻域中Log-Garbor作為基礎函數對圖像進行濾波,Log-garbor函數如下:

單演濾波方法是首先利用Log-Garbor對圖像進行濾波得到濾波后的信號f(x,y),之后利用濾波函數H對f(x,y)進行濾波得到濾波后的信號f h(x,y).將f(x,y)的實部fr(x,y),和f h(x,y)的實部f hr(x,y)和虛部f hi(x,y)組合到一起,我們就得到了一組二維的“解析信號”:

我們稱之為單演信號,在球坐標系下?∈[0,2π),θ∈[0,π),則有

因此單演信號的局部振頻(Local Amplitude)為:

局部相位(Local Phase)為:

局部方位(Local Orientation)為:

圖像總的振頻|E|、相位角?、方位角θ分別為:

這里n代表不同的尺度,這樣單演濾波的相位一致性就可以根據(2)式得到,圖像特征結構越明顯,它的相位一致性和總相位角也就越大,而方位角則分別反映圖像中這些特征結構的方向性.

圖2展示了一幅冕環(huán)圖像,該圖顯示的是搭載在“太陽動力學天文臺”(Solar Dynamics Observatory,SDO,2010)[32]上的“大氣成像組件”(Atmospheric Imaging Assembly, AIA)[33]在2015年6月18日00:00:11 UT 171?A波段觀測數據.該圖的空間分辨率為0.6′′,圖中數據已經經過了標準的AIA PREP校準過程.從圖2中我們可以發(fā)現(xiàn)多數冕環(huán)結構的亮度從足點到拱頂逐漸下降,導致冕環(huán)結構的清晰度從足點向上逐漸降低,最終在接近冕環(huán)頂部附近的冕環(huán)變得模糊難辨;另外我們也可以注意到大部分冕環(huán)結構有著一定的透明度,導致冕環(huán)的亮度變化更加復雜;還有部分冕環(huán)結構互相交叉,在它們交叉部分冕環(huán)的走向變得難以分辨,這都為冕環(huán)結構的識別帶來了巨大挑戰(zhàn).接下來我們分別利用高通濾波(high pass filtered)、巴特沃思帶通濾波(band pass filtered)、高斯帶通濾波、基于單演濾波的二維相位一致性方法對圖2中的測試圖像進行處理,結果如圖3所示.其中(a)、(b)、(c)分別是高通濾波、巴特沃思帶通濾波、高斯帶通濾波的結果,這里高通濾波是利用原圖減去它的16×16 boxcar平滑結果所得,巴特沃思帶通濾波、高斯帶通濾波使用IDL函數庫里面的BANDPASS FILTER程序實現(xiàn),其中巴特沃思帶通濾波的階數為20,巴特沃思和高斯帶通濾波的頻帶范圍為最高頻率的0.05–0.25倍.圖3的(d)、(e)、(f)為相位一致性方法所得到的相位一致性PC、相位角、方位角,計算過程中我們應用了單演濾波的方法,利用(8)式所給出的Log-Garbor函數作為基函數進行濾波,不同尺度濾波函數中的中心波長的取值為λ0=3×2.1n(這里中心波長的單位是像素),式中不同的n代表不同尺度的濾波函數,這里我們取n=0,1,2,3 4個尺度濾波函數.從結果我們可以看出相位一致性方法的確可以很好地對復雜環(huán)境下的冕環(huán)結構進行增強,尤其是相位角(如圖3(e)所示)可以很好地突出測試圖中的冕環(huán)結構.進一步對比表明在相位角圖像中整條冕環(huán)對比度明顯有著較大提高,尤其是冕環(huán)頂部較暗的部分對比度得到了很好的提升.

圖2 SDO衛(wèi)星171?A波段的冕環(huán)圖像.為了更加清晰地顯示圖像細節(jié),圖中顏色越黑代表輻射越強.虛線S0展示了圖4中的亮度變化曲線所在的位置.Fig.2 The coronal loop im ages observed by SDO 171?A.B lack m eans h igh fluxes in th is im age to show the details.T he dashed line S0 ind icates the location of the p rofiles in Fig.4.

為進一步對比不同圖像增強方法的效果,我們把幾種不同方法得到的S0(如圖2中黑色虛線所示)上的亮度變化曲線展示在圖4之中.如圖4所示S0中明顯有3個峰值分別代表3個冕環(huán)結構.從原圖亮度曲線(黑色實線)上我們可以看出,3個冕環(huán)中中間冕環(huán)有著較好的對比度,而兩邊的冕環(huán)則比較暗;從相位一致性相位角結果(綠色點線)我們可以看出,相位一致性方法對于S0兩邊的暗冕環(huán)結構的提取結果同中間亮冕環(huán)結構達到同一水平,而其他方法對暗結構的提取結果則明顯差于相位一致性的結果.總體而言,相位一致性相位角圖像不但可以很好地提取亮冕環(huán)結構,同時也可以很好地提取較暗冕環(huán)結構.最后我們注意到相位一致性方法也存在著一定的問題:它在增強了冕環(huán)結構的同時,也對其他的一些非冕環(huán)結構和噪聲進行了增強.

3 冕環(huán)識別算法

從我們的測試圖像可以看出,冕環(huán)結構有著亮度變化大,以及相互交疊等特點,這都對冕環(huán)識別提出了很大的挑戰(zhàn).但是通過大量的觀測分析我們也發(fā)現(xiàn)冕環(huán)結構基本上都是較光滑的大尺度結構,它們多數是從活動區(qū)磁極發(fā)出并且有著很好的連續(xù)性,盡管隨著高度增加亮度會有明顯降低,但是形態(tài)變化卻很平滑,主要由橢圓形的閉合冕環(huán)和放射狀的開放冕環(huán)構成,這些特征則對冕環(huán)的識別非常有利.

圖3 利用高通濾波、帶通濾波以及相位一致性算法進行冕環(huán)結構增強所得到的結果圖像Fig.3 T he coronal loop structu res strengthened by h igh pass filtered,band pass filtered,and phase congruency detective m ethods,resp ectively

為了更好地識別冕環(huán)結構,我們在ODM方法的基礎上發(fā)展了一個專門針對極紫外波段觀測到的冕環(huán)圖像進行優(yōu)化的算法.該算法主要是依據冕環(huán)形態(tài)變化比較平滑的特性對冕環(huán)識別過程進行優(yōu)化得到的.改進后的算法是通過冕環(huán)的方向性來指引識別過程,而冕環(huán)的方向則是通過計算以初始點為中心的上下兩個180°扇面上各方向亮度平均值的分布確定.這個扇面的半徑即我們的搜索半徑(簡稱rs),rs的取值應該大于冕環(huán)的寬度,依據不同的情況可取冕環(huán)寬度的幾倍,我們在計算過程中取35像素.基本的計算過程為,首先計算預期方向周圍180°扇面上各個方向的亮度平均值,將亮度大于最大亮度95%的方向進行平均從而得到該點的推進方向,并將識別結果沿該推進方向推進一段較小的距離(該距離就是推進步長,默認為2像素),之后循環(huán)這個過程從而對圖像中的冕環(huán)結構逐步識別.直接使用這個算法則會遇到與ODM相同的問題,使得識別過程很容易被噪聲打斷,從而導致識別的結構不完整.針對這個問題,我們對上述算法進行了一些改進,其詳細過程如圖5所示.圖中兩個交叉的黑色弧形環(huán)代表我們要識別的冕環(huán)結構,其中一條環(huán)上的白色星線為我們已經識別的部分冕環(huán)結構,已經識別的冕環(huán)的切線方向為該點的預期方向(如圖5中綠色箭頭所示).這樣如果計算過程進行到冕環(huán)交叉的地方,或者冕環(huán)外面有一些亮度較大的噪聲點時,我們的算法將會識別出兩處或多處亮度較大的方向,圖中紅色扇面就標出了兩處不連續(xù)的亮度較大方向.如果將這兩處方向共同進行平均就很容易得到一個錯誤的方向,使推進點脫離冕環(huán)結構,最終導致識別提前結束.這種情況下由于冕環(huán)的形態(tài)變化比較平緩,因此下一點上冕環(huán)方向不會有太大的變化,所以這里我們只取預期方向周圍較小范圍(計算過程中默認為±10°)內的亮度滿足條件的方向進行平均.如果我們預期的結果可能出現(xiàn)的范圍內沒有滿足條件的方向或者整個搜索范圍內找到的亮度滿足大于最大亮度95%的方向總的個數超過一定范圍(我們計算過程中默認為方向總個數的1/3),則按照預期方向進行推進并且將該次識別記為一次異常識別.

圖4 S0上的亮度變化曲線,其中黑色實線為原始數據,藍、綠、紅色3條點線分別為利用單演濾波得到的相位一致性PC、相位角圖像以及高通濾波所得圖像的結果,而黑色點線以及虛線則是巴特沃思帶通濾波和高斯帶通濾波所得到圖像的結果,這里所有結果都被等比例歸一化到0–255.F ig.4 The in tensity p rofiles of S0.The b lack so lid line ind icates the origina l data.T he b lue,green,and red dotted lines show the data got from PC im age,im age,and high pass filtered im age,respectively. The Bu tterw orth and G ausssian band pass filtered resu lts are show n by the b lack dotted line and the b lack dashed line,resp ectively.Notew orthily,all the resu lts are norm alized to 0–255.

圖5 冕環(huán)識別算法示意圖.其中綠色箭頭為預期的冕環(huán)方向,其周圍較小的綠色扇面為我們程序所允許的結果方向的范圍,藍色180°扇面則是我們的搜索范圍,紅色扇面則是具有較大平均亮度的兩個方向范圍.Fig.5 The schem atic of ou r loop-tracing m ethod.The green arrow ind icates the expected d irection. A round th is green arrow,the sm aller green sector and 180°b lue sector are the range of p ossib le resu lt a llow ed by ou r p rogram and the search ing scop e,resp ectively.T he red sectors ind icate the tw o ranges of d irection along w h ich p ixels have a larger m ean intensity.

4 冕環(huán)自動識別方法

基于相位一致性所得到的冕環(huán)增強圖像和經過改進的冕環(huán)識別算法,我們提出了一個新的冕環(huán)自動識別方法(Phase-congruency-based method,簡稱PCB方法).這個方法主要通過兩個步驟進行冕環(huán)自動識別.首先我們通過一個程序對原始日冕觀測圖像進行預處理,該過程要先對原始圖像進行分塊,每一塊子圖的平均值要大于一定的閾值mi(我們的計算中默認取值為0.05,代表圖像最大值的0.05倍),則認為該塊圖像中可能存在冕環(huán)結構.如果存在冕環(huán)結構則將該塊子圖中的最大值所在的位置標記為一個初始點,這些點將被用作下一步冕環(huán)結構識別的初始點,接下來該程序則會應用相位一致性方法對原始圖像進行冕環(huán)結構增強.經過實驗,我們發(fā)現(xiàn)針對AIA和TRACE望遠鏡171?A波段的觀測數據,取λ0=3×2.1n(n=0,1,2,3)像素的4個不同尺度的濾波函數進行計算所得到的相位角圖像,即可使絕大部分的冕環(huán)結構得到很好的增強,因此我們將以該圖為參考進行第2步的冕環(huán)識別.第2步,我們將依據第1步所得到的冕環(huán)增強圖像和初始點進行冕環(huán)識別,識別過程中主要用到我們改進過的冕環(huán)識別算法對每一個初始點進行冕環(huán)的識別.這里我們首次將識別方向的改變情況以及識別結果點的亮度共同用于判斷識別結果的正確性,如果連續(xù)zn次識別結果(默認zn=2)均為異常識別或者識別結果點的亮度低于我們提前設定的截止閾值toff(默認為0.3,代表圖像最大亮度的0.3倍),則認為識別已經超出真正的冕環(huán)結構并且終止本次識別.最后我們的程序將識別的長度大于30像素的結構進行記錄,而更短的結構則認為是其他結構,不進行任何記錄.

另外,我們也注意到前面的程序將可能在同一個冕環(huán)結構上產生多個初始點,這樣將會導致同一冕環(huán)被重復識別多次,在我們識別方法中提供了兩個辦法解決這個問題.第1種方法(簡稱為PCB1)是每次識別開始時我們就判斷該初始點與前面所識別出來的冕環(huán)上所有點的距離的最小值,如果小于我們設定的閾值(簡稱mw,默認為25像素),我們就將該初始點剔除,不再進行識別.考慮到第1種方法中如果閾值mw設定得太小,則還是會有少量重復識別結果,設定得太大時,可能錯誤地剔除部分靠得較近的相平行的冕環(huán),因此我們提出了第2種方法(簡稱PCB2).該方法是將每一個新識別的點同前面所有結果的距離的最小值都計算出來,如果小于我們設定的另一個閾值(簡稱md,默認為6像素)的點數超過3個則終止識別.這里以圖3中的測試圖像為例,以默認的參數進行識別,第1種方法中的mw取值為25,第2種方法中的md則取為6,它們的結果如圖6所示.圖中星狀點代表我們所得到的初始點,粗線代表程序識別的冕環(huán)結構,兩種方法在相同電腦上所用時間分別為3.463 s和8.502 s.經過對比我們可以發(fā)現(xiàn)這兩種方法基本上都可以有效地去除重復識別的結果,但是從圖6中上圖的效果看第1種方法的結果中還有少量重復識別發(fā)生,因此剔除效果可能稍微差于第2種方法,但是從計算速度上看其比第2種方法要快.

值得注意的是,終止識別的條件的選取一直是一個復雜的問題,如果取得過于苛刻則不能很好地識別較暗的冕環(huán)和一些冕環(huán)中較暗的部分,如果取得過于寬松則亮結構冕環(huán)的識別容易超出其邊界,得到錯誤的識別結果,這導致準確性和完整性很難同時都獲得.但是由于超出冕環(huán)后的錯誤識別結果不再具有很好的連續(xù)性,所以超出后的識別方向會有劇烈的波動.根據這些特點我們首次將識別方向的變化情況用于判斷是否停止識別,這樣我們就可以進一步將原來的條件放寬松.試驗表明較小的亮度截止閾值toff加上較小的方向改變閾值zn可以很好地彌補之前較寬松的截止亮度閾值產生的錯誤識別結果,而且可以得到更加完整平滑的冕環(huán)結構,這些都使得我們的識別結果更加接近真實情況.

5 識別結果的比較和分析

針對冕環(huán)的特性和相位一致性方法所得到的冕環(huán)增強圖像的特點,我們提出了一個新的PCB冕環(huán)識別方法.為了更好地同其他方法進行對比,這里我們采用Aschwanden等[21]和M cAteer等[25]所用的測試數據對我們的方法進行測試,測試中所用的參數均采用前文給出的默認值.TRACE的這組數據如圖7所示,該數據是TRACE衛(wèi)星在1998年5月19日22:21 UT 171?A波段的觀測結果.利用PCB2方法所得到的識別結果展示在圖8之中,圖中我們可以看到PCB2識別的冕環(huán)整體上都比較完整,而且識別的冕環(huán)結構都比較光滑,尤其是對圖像兩邊的開放冕環(huán)的識別有著很大的優(yōu)勢,接下來我們通過一些定量參數以及局部區(qū)域結果對PCB2的具體特性進行詳細分析.

圖6 兩種不同的避免重復方法的識別結果.上下兩圖分別顯示PCB 1和PCB 2所識別的結果.Fig.6 T he iden tified resu lts of tw o m ethods that are used to avoid rep eating iden tifications.Top (bottom)panel is the resu lt of PCB 1(PCB 2).

5.1 完整性分析

冕環(huán)識別結果的完整性可以分為兩個方面:首先是識別的冕環(huán)種類完整,其次就是識別出來的冕環(huán)的長度完整.由于各個不同方法的識別結果都含有大量的冕環(huán),而且識別結果的位置也各有不同,所以很難對其完整性進行直接對比.根據Aschwanden等[21]所采用的一種方法,我們可以通過識別的長度大于L的冕環(huán)總數隨L的分布(簡稱冕環(huán)累積分布)來大致反映冕環(huán)識別結果的完整度.該參數既反映了識別冕環(huán)的長度,又反映了冕環(huán)的數量.為了對比PCB2方法識別結果的完整性,我們也計算出了我們的PCB2方法識別結果的冕環(huán)累積分布,分布結果展示在圖9之中.圖像表明我們的方法識別出了96條長度大于70像素的冕環(huán),與Aschwanden等[21]探測到大于70像素冕環(huán)最多的RAS方法所得到的結果91條水平相當.除此之外,為了更全面地反映我們的方法識別冕環(huán)的完整性,我們也將識別結果同Aschwanden等[21]的結果(圖3、圖4、圖5)進行詳細對比.一方面從識別的種類上看,PCB2方法所識別的圖像中間部分的閉合冕環(huán)和兩側的開放冕環(huán)幾乎同Aschwanden等[21]所提到的OCM、DAM、UDM、ODM 4種方法的總和相當;另一方面在識別冕環(huán)的長度上,我們的識別結果在閉合冕環(huán)的頂部并沒有像UDM和ODM那樣發(fā)生中斷,而在兩邊開放冕環(huán)的頂端,我們的方法則可以將大多數開放冕環(huán)識別到接近圖像邊緣,這方面PCB2方法對開放冕環(huán)的探測結果甚至比Aschwanden等[21]所提到的5種方法都要好.

圖7 TRACE 171?A波段于1998年5月19日22:21 UT觀測數據.與前面A IA的圖像相同,圖中顏色越黑代表輻射越強.Fig.7 T he corona l im age observed by the TRACE 171?A at 22:21 UT on 1998 M ay 19.Notew orthily, b lack m eans h igh fluxes in th is im age.

圖8 TRACE數據的冕環(huán)自動識別結果Fig.8 The coronal loops autom atically ex tracted from the tested TRACE im ages

圖9 TRACE數據識別中冕環(huán)長度大于L的累計分布圖Fig.9 T he cum u lative d istribu tion of loops iden tified from the tested TRACE im age with a length larger than L

5.2 準確性和靈敏度分析

為了對冕環(huán)識別結果的準確性和靈敏度進行詳細的比較分析,我們選取了圖像中A、B兩個小區(qū)域來進行詳細對比,該區(qū)域同Aschwanden等[21]所選取的區(qū)域完全相同.圖10和圖11分別展示了我們選取的A、B兩個區(qū)域,兩幅圖的背景是該區(qū)域高通濾波圖像的等高線,里面的粗線即為PCB2方法的識別結果,圖中等高線密集的地方就是冕環(huán)所在的位置.通過對比識別結果和等高線的位置關系,我們就可以知道識別結果是否準確,如果識別的冕環(huán)和背景等高線分布密集的地方符合得較好,則說明識別的結果比較準確.從圖10和圖11的結果我們可以看到,PCB2方法識別冕環(huán)結構與背景等高線很好地符合在一起.同文獻[21]中的結果相比,可以看出我們的結果有以下特點:首先,PCB2的結果同OCM和DAM比較類似,比較光滑,不像UDM那樣識別結果有著許多小的虛假振蕩;其次,結果冕環(huán)末端也比較平滑,這一點要優(yōu)于RAS的結果;最后,識別結果中錯誤識別水平也很低.通過對圖像中子區(qū)域識別結果仔細對比分析,我們可以看出PCB2冕環(huán)識別方法的結果同OCM、DAM、UDM、ODM、RAS相比處于一個較準確的水平.而靈敏度方面,PCB2方法可以識別出兩個區(qū)域中絕大多數的冕環(huán)結構,識別的數目和完整性同靈敏度較高的RAS處于同一水平,但是準確性方面卻彌補了RAS識別的冕環(huán)末端存在少量錯誤的彎曲和錯誤識別一些小尺度的非冕環(huán)結構兩方面的不足.由于PCB2方法在完整性以及準確性方面的優(yōu)勢,尤其是對開放冕環(huán)識別方面的良好結果,我們可以發(fā)現(xiàn)PCB2方法的確是一個有效可行的冕環(huán)識別方法.

圖10 TRACE數據A區(qū)域的詳細識別結果Fig.10 The detailed ex traction resu lt of the TRACE tested data in region A

圖11 TRACE數據B區(qū)域的詳細識別結果Fig.11 The detailed extraction resu lt of the TRACE tested data in region B

6 典型樣例中的應用和分析

由于SDO衛(wèi)星的發(fā)射,使得我們可以獲得大量的日冕觀測數據,這里我們也將PCB2方法用于更多的SDO/AIA 171?A波段的觀測數據之中.首先我們選取了4個典型的冕環(huán)觀測數據用來測試和分析我們的方法,數據的詳細信息如表1所示,測試之前幾組數據已經經過標準的預處理過程.數據1–4在撒點過程中的閾值mi分別取為0.03、0.02、0.03、0.02(該閾值的選取既要保證識別較多的結果,又不能包含太多的虛假結構),冕環(huán)識別過程中所需要的搜索半徑rs和推進步長,以及防止重復識別的閾值md,還有用于停止識別的閾值toff和zn分別取前面給出的默認值.對表1中數據進行冕環(huán)識別的結果如圖12和圖13所示,其中圖12上下兩圖分別展示了測試數據1和數據2的結果,測試數據3和測試數據4的結果則分別展示在圖13中的上下兩圖.在每一幅結果圖像中,粗線就是PCB2方法的識別結果,圖像的背景為原始觀測數據,疊加在背景上的等高線與圖10和圖11相同,為高通濾波后圖像的等高線.與本文第5章的方法相同,我們可以通過分析圖像中閉合冕環(huán)和開放冕環(huán)頂端的識別情況來判斷識別結果是否完整,對比結果圖像中等高線的位置和識別冕環(huán)的對應情況來反映識別結果的準確性.首先我們可以發(fā)現(xiàn)PCB2方法的識別結果的確和圖像等高線的分布有著很好的對應,因此可以確定PCB2方法的識別結果有著很好的準確性,其次從識別種類上看,不管是較直的開放冕環(huán),還是有一定弧度的閉合冕環(huán),都可以被PCB2方法所識別.通過對識別出的冕環(huán)結果進行分析,我們還可以看到不管是較亮的足點部分還是較暗的冕環(huán)頂部都可以被較好地識別.但是我們也注意到結果中存在少量將其他結構誤識別為冕環(huán)的情況,這些誤識別結構是位于冕環(huán)外圍譜斑附近一些大尺度并且較狹窄的亮區(qū)域,這些區(qū)域同較暗的冕環(huán)結構很相似,因此比較容易被誤識別.

表1 測試圖像的詳細信息Tab le 1 The detailed in form ation of the test im ages

7 總結與討論

冕環(huán)是日冕中普遍存在的由磁場束縛的熱等離子組成的亮環(huán)狀結構[1],它很好地反映了我們平時很難直接觀測到的日冕中的磁場結構[34].大量的觀測數據表明多數冕環(huán)都處在一個復雜的環(huán)境中,而且其本身的亮度也有著很大變化,這都為我們識別冕環(huán)結構帶來了很大的挑戰(zhàn).針對冕環(huán)結構的自動識別,目前已經有一些辦法被提出,這些方法都可以在一定程度上得到較好的結果,但是也存在著一些不足.這些冕環(huán)識別的方法一般包含對原圖像中的冕環(huán)結構進行增強,以及對增強結果進行冕環(huán)標記識別兩個過程.而在2014年Feng等[26]就對利用相位一致性方法來增強日冕中冕環(huán)結構的方法進行了詳細研究,結果表明相位一致性的確對日冕中的冕環(huán)結構有著很好的增強作用,但是基于相位一致性冕環(huán)增強圖像的冕環(huán)識別方法卻一直沒有被提出.本文我們將相位一致性中單演濾波方法所得到冕環(huán)增強圖像用于冕環(huán)識別之中,提出了一種基于相位一致性的冕環(huán)自動識別方法.同時我們也將冕環(huán)結構識別算法依據冕環(huán)形態(tài)變化比較平滑的結構特性進行了改進,只取預期方向周圍較小范圍內的方向作為有效方向來計算真正的推進方向.最后我們對多組冕環(huán)的觀測數據進行測試和分析,測試的結果表明我們的方法的識別結果不但具有很高的完整性,而且同時具有較高的準確性,這是大多數方法所不能同時具備的.而且測試結果也表明PCB2方法對大尺度冕環(huán)結構有著很好的識別效果,尤其在大尺度開放冕環(huán)結構的識別方面有著明顯的優(yōu)勢.以上結果表明相位一致性方法增強的冕環(huán)圖像的確可以被用于冕環(huán)結構識別,而我們針對相位一致性增強圖像所開發(fā)的冕環(huán)自動識別方法也確實可以得到比較好的識別結果,是一套現(xiàn)實可行的冕環(huán)自動識別方法,對冕環(huán)的探索和研究都有一定的使用價值.

圖12 PCB 2方法在A IA數據樣本1(上)和2(下)中識別的冕環(huán)結果Fig.12 The loops detected by PCB 2 from the A IA sam p les 1(top)and 2(bottom)

圖13 PCB 2方法在A IA數據樣本3(上)和4(下)中識別的冕環(huán)結果Fig.13 The loops detected by PCB 2 from the A IA sam p les 3(top)and 4(bottom)

此外我們也注意到,相位一致性方法增強的日冕觀測圖像中的一些噪聲和非冕環(huán)結構也同時被增強,針對這樣的情況我們提出了一個綜合的冕環(huán)自動識別方法來減少這些被增強的噪聲和非冕環(huán)結構對識別過程的影響.在PCB算法中我們主要通過3個方面共同作用來避免識別圖像中的非冕環(huán)結構.首先我們選取初始點時就對初始點周圍一個小范圍內的平均值有一定的要求,平均值要大于我們設定的閾值mi才有可能成為一個初始點,這是因為冕環(huán)觀測圖像中真正的冕環(huán)結構通常都比較亮,而其他一些亮斑或者非冕環(huán)結構周圍都比較暗,所以閾值mi的設定可以幫助我們有效地排除多數的噪聲以及非冕環(huán)結構.第2個方面就是我們將識別方向的改變情況以及識別結果點的亮度共同用做停止識別的判據,如果連續(xù)zn次方向識別結果均為異?；蛘咦R別點的亮度小于閾值toff則終止識別.這主要是相位一致性增強后的冕環(huán)結構有著很高的對比度和連續(xù)性,而其他的噪聲或者非冕環(huán)結構則不同,相比冕環(huán)結構圖像中的冕環(huán)結構,它們一般沒有規(guī)則的形狀,連續(xù)性較差,方向改變也較大.因此zn和toff兩個判據的共同限制也可以幫助我們排除大部分的非冕環(huán)結構.最后我們還要將識別結果的長度做為一個標準來降低錯誤識別的可能性,這是因為真正的冕環(huán)結構屬于較大尺度的結構,因此在長度上要遠長于大部分的其他結構.值得注意的是,雖然經過這3個步驟的共同作用已經將錯誤識別的概率降到了很低的水平,但是依然可能有極少量的非冕環(huán)結構被識別,這些結構同活動區(qū)外邊緣較暗的冕環(huán)結構很相似,因此很難被區(qū)分.

最后應該指出的是,雖然我們提出的冕環(huán)識別算法主要是針對相位一致性方法進行優(yōu)化,但是這個算法也為其他方法的優(yōu)化提供了參考.我們的結果表明,冕環(huán)結構本身的特性為我們進一步優(yōu)化冕環(huán)識別方法提供了一個很好的參考.而且我們也注意到我們的方法也存在少量錯誤識別等一些問題,因此在這些方面還需要更深入的研究.

[1]Vaiana G S,Reidy W P,Zehnp fennig T,et al.Science,1968,161:564

[2]M oortel I,Nakariakov V M.RSPTA,2012,370:3193

[3]Parnell C E,De M oortel I.RSPTA,2012,370:3217

[4]Go lub L,Noci G,Poletto G,et a l.A p J,1982,259:359

[5]Vaiana G S,K rieger A S,T im othy A F.SoPh,1973,32:81

[6]Foukal P V.A p J,1976,210:575

[7]B rekke P,K jeldseth-M oe O,Harrison R A.SoPh,1997,175:511

[8]Huang Z,Frazin R A,Land i E,et a l.A p J,2012,755:86

[9]Rea le F.LRSP,2014,11:4

[10]A schw anden M J,N itta N V,W¨u lser J P.A p J,2008,680:1477

[11]Go lub L,Nystrom G,Heran t M,et al.Natu re,1990,344:842

[12]Kaiser M L,Kucera T A,Davila J M,et al.SSRv,2008,136:5

[13]張雪飛,劉煜,申遠燈,等.天文學進展,2012,30:159

[14]U ritsky V M,Dav ila J M,V iall N M,et a l.A p J,2013,778:26

[15]李海東,趙麗,梁紅飛,等.天文學報,2013,54:27

[16]Li H D,Zhao L,Liang H F,et a l.ChA&A,2013,37:266

[17]W ang T,ofm an L,Davila J M.Ap J,2013,775:L23

[18]張真,何馨,梁紅飛.天文學報,2013,54:39

[19]Zhang Z,He X,Liang H F.ChA&A,2013,37:428

[20]Lee J K,Newm an T S,G ary G A.Pattern Recogn it,2006,39:246

[21]A schw anden M J,Lee J K,G ary G A,et a l.SoPh,2008,248:359

[22]Inhester B,Feng L,W iegelm ann T.SoPh,2008,248:379

[23]Parent P,Zucker S.IEEE Trans,1989,11:823

[24]Handy B N,A cton L W,K ankelborg C C,et a l.SoPh,1999,187:229

[25]M cA teer R T J,K estener P,A rneodo A,et a l.SoPh,2010,262:387

[26]Feng S,Xu Z,W ang F,et a l.SoPh,2014,289:3985

[27]Oppenheim A V,Lim J S.P roceed ings of the IEEE,1981,69:529

[28]M orrone M C,Ow ens R.Pattern Recogn ition Lett.,1987,6:303

[29]Kovesi P.V idere:Journal of Com pu ter V ision Research,1999,1:1

[30]K ovesi P.Psycho logica l Research,2000,64:136

[31]Felsberg M,Somm er G.M ustererkennung 2000,22.DAGM-Sym p osium.London:Sp ringer-Verlag, 2000:195-202

[32]Pesnell W D,T hom pson B J,Cham berlin P C.SoPh,2012,275:3

[33]Lem en J R,T itle A M,A k in D J,et a l.SoPh,2012,275:17

[34]Rosner R,Tucker W H,Vaiana G S.A p J,1978,220:643

A N ew Coronal Loop Identified Technology Based on the Phase Congruency Detective M ethod

LIHong-bo1,2ZHAO M ing-yu1LIU Yu1

(1 Yunnan O bserva to ries,Chinese A cadem y of Scien ces,K unm ing 650011) (2 Un iversity of Chinese A cadem y of Scien ces,Beijing 100049)

We test a phase congruency method for automatically extracting coronal loops from the coronal observation data,and propose a new coronal loop detective method.On account of smooth changes of the direction along coronal loops,the possib le direction of loops is restricted in a sm all range for im proving the detective result. Beyond that,inspired by the structural characteristics of coronal loops,we firstly suggest that both the change of result directions and the flux of resu lt points can be used to term inate the loop detection.Finally,several coronal images are detected by our loop extracting method,and the result indicates that the congruency enhanced im age is exactly suitable for coronal detection and extraction,and our new loop detective method can also extract most of loops in the test im ages com p letely and accurately. Besides that,our research also suggests that the characteristics of coronal loop should be a significant reference for im p roving other loop detectivemethods.

sun:corona,techniques:im age processing,methods:observational

P182;

A

10.15940/j.cnki.0001-5245.2016.04.003

2015-11-24收到原稿,2016-03-15收到修改稿

?國家自然科學基金項目(11533009)資助

?hbli@ynao.ac.cn

?m yzhao@ynao.ac.cn