太陽射電II型爆發(fā)中頻率分裂現(xiàn)象的特征自動(dòng)提取研究

李國良，袁國武，高冠男

(1.中國林業(yè)科學(xué)研究院資源昆蟲研究所,云南昆明，650233;2.云南大學(xué)信息學(xué)院,云南昆明，650091;3.中國科學(xué)院云南天文臺,云南昆明，650011)

0 引言

隨著太陽觀測儀器的快速發(fā)展，大量的觀測數(shù)據(jù)被采集，但是反映太陽爆發(fā)活動(dòng)的數(shù)據(jù)只占所采集數(shù)據(jù)的1%[1]。對于觀測數(shù)據(jù)，天文學(xué)家通常對太陽射電爆發(fā)的部分更感興趣，而在海量的觀測數(shù)據(jù)中，找出太陽射電爆發(fā)就成了一項(xiàng)繁重、枯燥的工作，如果再進(jìn)一步對爆發(fā)中包含的精細(xì)結(jié)構(gòu)進(jìn)行人工提取，就更加重了天文學(xué)家的負(fù)擔(dān)。因此，對太陽射電頻譜儀采集的數(shù)據(jù)進(jìn)行自動(dòng)射電爆發(fā)檢測和精細(xì)結(jié)構(gòu)的特征自動(dòng)提取，對于太陽射電天文研究非常有意義。

太陽射電頻譜儀對太陽射電進(jìn)行觀測，得到的數(shù)據(jù)是強(qiáng)度隨輻射頻率和時(shí)間變化的圖譜，稱為射電動(dòng)態(tài)頻譜圖。根據(jù)動(dòng)態(tài)頻譜上的射電信號變化，我們可以找到頻率隨時(shí)間緩慢變化的II型太陽射電暴。1947年，Payne-Scott[2]等人在200MHz、100MHz和60MHz頻率上發(fā)現(xiàn)了這種具有慢速頻率漂移的太陽射電爆發(fā)。Wild和McCready[3]于1950年將其定義為II型射電暴。經(jīng)過多年觀測發(fā)現(xiàn)II型射電暴通常都展現(xiàn)出基頻和諧頻輻射結(jié)構(gòu)[4]，在20%的II型爆發(fā)中還發(fā)現(xiàn)基頻或二次諧頻會(huì)分裂成形態(tài)相似的兩個(gè)部分這一頻率分裂現(xiàn)象[5]，該現(xiàn)象的特征參數(shù)對于理解II型射電暴的產(chǎn)生機(jī)制十分重要。近十幾年來，國內(nèi)外學(xué)者僅對少數(shù)精細(xì)結(jié)構(gòu)開展了特征提取研究，尚未涉及頻率分裂現(xiàn)象。所以，目前開展頻率分裂現(xiàn)象的特征自動(dòng)提取研究將會(huì)對后續(xù)的太陽物理研究工作提供很大幫助。

1 本文采用的太陽射電頻譜數(shù)據(jù)

中國科學(xué)院云南天文臺擁有兩臺具有了極高頻率分辨率和時(shí)間分辨率的全數(shù)字化太陽射電頻譜儀，能夠?qū)Ω鼮榫?xì)和快速變化的射電爆發(fā)及精細(xì)結(jié)構(gòu)進(jìn)行觀測。本文實(shí)驗(yàn)所采用的太陽射電頻譜數(shù)據(jù)來自云南天文臺撫仙湖太陽觀測基地的11米口徑米波太陽射電頻譜儀，這臺米波射電頻譜儀于2011年7月開始運(yùn)行，工作頻段為70-700MHz，最高頻率分辨率為200kHz，優(yōu)于德國AIP的幾百KHz和澳大利亞Culgoora的1MHz分辨率[6]，處于國際領(lǐng)先水平?？紤]到時(shí)間分辨率和頻率分辨率之間的矛盾，該儀器日常工作的時(shí)間分辨率被設(shè)置為80ms，頻率分辨率設(shè)置為200kHz。工作時(shí)間為北京時(shí)間8:00-18:00。

表1 云南天文臺米波太陽射電頻譜儀的觀測參數(shù)

2 II型射電爆發(fā)區(qū)域的確定

2.1 太陽射電頻譜圖像去噪

本文對原始頻譜中的通道做歸一化處理，消除強(qiáng)背景干擾和去除通道差異產(chǎn)生的噪聲[7]，增強(qiáng)圖像中有用的信息。通道歸一化的原理是用頻譜圖像中每個(gè)通道的凈太陽射電爆發(fā)流量來除以該通道的平均值，以此減小各個(gè)通道之間的差異來實(shí)現(xiàn)噪聲消除。由于單通道的凈太陽射電爆發(fā)流量?fy(x)為原始通道的流量fy(x)減去該通道的平均值，所以單通道的歸一化公式為通過對y取全通道信息實(shí)現(xiàn)整幅圖像的通道歸一化。

2.2 大津法二值化

大津法即最大類間方差法，能夠根據(jù)太陽射電頻譜圖像本身的灰度特性,將圖像分割為目標(biāo)前景與背景兩個(gè)部分。算法公式為分別為前景點(diǎn)數(shù)和背景點(diǎn)數(shù)所占圖像的比例，μa和μb分別為前景和背景的平均灰度。算法能自動(dòng)選取最佳閾值T計(jì)算目標(biāo)前景與背景的最大類間總方差(T)，實(shí)現(xiàn)二值分割。

2.3 形態(tài)學(xué)重建確定爆發(fā)區(qū)域

3 頻率分裂現(xiàn)象的特征自動(dòng)提取

經(jīng)典的精細(xì)結(jié)構(gòu)特征提取方法有連片搜索算法[9]、level set提取法[10]、小波變換分離法[11]，這些方法僅適用于區(qū)域型精細(xì)結(jié)構(gòu)的提取。針對頻率分裂現(xiàn)象這類邊緣型的精細(xì)結(jié)構(gòu)，本文提出了新的方法。

圖1 爆發(fā)區(qū)域確定實(shí)驗(yàn)

3.1 高斯卷積平滑圖像

高斯濾波是一種線性平滑濾波，常用于圖像的平滑處理。本文使用一個(gè)標(biāo)準(zhǔn)差為0.5，大小為5×5的卷積核對爆發(fā)區(qū)域進(jìn)行二維高斯卷積平滑處理。對頻譜圖像進(jìn)行平滑處理可以降低爆發(fā)區(qū)域噪聲，使爆發(fā)邊界更加清晰，增強(qiáng)精細(xì)結(jié)構(gòu)的特征。

3.2 爆發(fā)區(qū)域的邊緣檢測

為了消除不良干擾,保留爆發(fā)區(qū)域的重要結(jié)構(gòu)信息，需要對平滑處理后的爆發(fā)進(jìn)行邊緣檢測。文獻(xiàn)[12]中，P.J.Zhang等人利用Canny邊緣檢測算法對來自NDA的射電數(shù)據(jù)進(jìn)行邊緣檢測，得到很好的實(shí)驗(yàn)效果。但是采用Canny算法對云南天文臺的數(shù)據(jù)進(jìn)行處理會(huì)加強(qiáng)噪聲，而且檢測出的邊緣呈不閉合狀態(tài)，其他傳統(tǒng)邊緣檢測方法的效果同樣不佳，因此本文根據(jù)邊緣指示函數(shù)[13]設(shè)計(jì)了適合的算法。

3.3 基于距離正則化水平集方法自動(dòng)提取頻率分裂現(xiàn)象

距離正則化水平集是以梯度流形式進(jìn)行演化的，能夠有效處理分裂與合并這一類型的拓?fù)渥兓痆14]，本文基于Li[14]等人提出的DRLSE算法進(jìn)行編程實(shí)現(xiàn)并應(yīng)用在基于邊緣的主動(dòng)輪廓模型上，實(shí)現(xiàn)距離正則化水平集方法自動(dòng)提取精細(xì)結(jié)構(gòu)。針對頻率分裂現(xiàn)象的不規(guī)則性，該方法能很好的實(shí)現(xiàn)輪廓特征提取，而且演化過程非常穩(wěn)定。

（1）初始輪廓設(shè)置

太陽射電II型爆發(fā)發(fā)生在頻率較低的區(qū)域，起始頻率不會(huì)超過幾百兆赫茲。2007年White S針對觀測到的頻率分裂高頻事件報(bào)道了頻率分裂現(xiàn)象的基頻起始頻率為380MHz，二次諧頻的起始頻率為760MHz[6]這一研究成果，這為頻率分裂現(xiàn)象的觀測和進(jìn)一步研究提供了一大依據(jù)。因此，結(jié)合本文實(shí)驗(yàn)數(shù)據(jù)自身的頻段信息，將初始輪廓位置設(shè)置在380MHz到700MHz的范圍內(nèi)。

（2）頻率分裂現(xiàn)象的自動(dòng)提取

距離正則化水平集方法依賴于圖像的梯度和能量函數(shù)實(shí)現(xiàn)演化功能，能量函數(shù)最小化運(yùn)算驅(qū)動(dòng)輪廓以梯度流形式進(jìn)行演化，圖2（c）中的等高線即為能量值變化的體現(xiàn)，初始輪廓會(huì)向著梯度較大的邊緣不斷進(jìn)行能量最小化運(yùn)算以接近邊緣。當(dāng)輪廓越來越接近邊緣時(shí)，算法會(huì)減緩輪廓的演化直至停止。本文根據(jù)經(jīng)驗(yàn)將演化的迭代次數(shù)設(shè)為固定值1200，能夠滿足云南天文臺數(shù)據(jù)的精細(xì)結(jié)構(gòu)提取。

實(shí) 驗(yàn)定義能量函數(shù)為E(?) =μRp(?) +Eext(?)，其中μRp(?)為正則項(xiàng)，能夠保持演化的穩(wěn)定性。外部能量函數(shù)Eext(?)中包含邊緣項(xiàng)和區(qū)域向 ：Eext(φ) =λLg(φ) +αA(φ)，通過最小化外部能量函數(shù)能夠?qū)崿F(xiàn)演化功能。其中，區(qū)域項(xiàng)αA(φ)能夠確定演化方向，邊緣項(xiàng)λLg(φ)取最小值時(shí)說明輪廓已經(jīng)演化到目標(biāo)邊界處，演化停止。

（3）精細(xì)結(jié)構(gòu)提取結(jié)果及分析

實(shí)驗(yàn)得到如圖2（d）中所示的提取結(jié)果，說明本文方法確實(shí)能有效實(shí)現(xiàn)頻率分裂現(xiàn)象的特征自動(dòng)提取。由于本文數(shù)據(jù)中存在破壞了II型射電暴完整性的壞通道，導(dǎo)致輪廓在垂直方向上的演化最終停在射電暴的截?cái)嗵?，但這并不影響最終的提取效果。假設(shè)沒有壞通道的干擾，輪廓在迭代1200次后同樣會(huì)在恰當(dāng)?shù)奈恢猛Ｖ?，最終提取出頻率分裂現(xiàn)象，而不會(huì)肆意演化。

圖2 頻率分裂現(xiàn)象提取實(shí)驗(yàn)

精細(xì)結(jié)構(gòu)中的大部分特征參數(shù)只與其所在位置有關(guān)與爆發(fā)強(qiáng)度無關(guān)[15]。所以在提取出頻率分裂現(xiàn)象后，可以通過記錄頻率分裂現(xiàn)象的頻率最大值點(diǎn)和最小值點(diǎn)所在的通道和時(shí)間來計(jì)算精細(xì)結(jié)構(gòu)的持續(xù)時(shí)間和帶寬，然后利用線性回歸方法計(jì)算出頻率漂移率[16]。

4 總結(jié)與展望

本文的不足之處在于滿足研究要求的II型爆發(fā)數(shù)據(jù)不多，所以在提取出精細(xì)結(jié)構(gòu)之后只是簡單描述了可行的特征參數(shù)統(tǒng)計(jì)分析方法，希望以后在更多數(shù)據(jù)積累下，實(shí)現(xiàn)完善的特征參數(shù)統(tǒng)計(jì)分析工作。另外，II型射電爆發(fā)中還存在Herringbone精細(xì)結(jié)構(gòu)，常常出現(xiàn)在頻率分裂現(xiàn)象的基頻和二次諧頻上，在本文的基礎(chǔ)上，可以開展相關(guān)研究。