李國良,袁國武,高冠男
(1.中國林業(yè)科學(xué)研究院資源昆蟲研究所,云南昆明,650233;2.云南大學(xué)信息學(xué)院,云南昆明,650091;3.中國科學(xué)院云南天文臺,云南昆明,650011)
隨著太陽觀測儀器的快速發(fā)展,大量的觀測數(shù)據(jù)被采集,但是反映太陽爆發(fā)活動的數(shù)據(jù)只占所采集數(shù)據(jù)的1%[1]。對于觀測數(shù)據(jù),天文學(xué)家通常對太陽射電爆發(fā)的部分更感興趣,而在海量的觀測數(shù)據(jù)中,找出太陽射電爆發(fā)就成了一項繁重、枯燥的工作,如果再進一步對爆發(fā)中包含的精細結(jié)構(gòu)進行人工提取,就更加重了天文學(xué)家的負(fù)擔(dān)。因此,對太陽射電頻譜儀采集的數(shù)據(jù)進行自動射電爆發(fā)檢測和精細結(jié)構(gòu)的特征自動提取,對于太陽射電天文研究非常有意義。
太陽射電頻譜儀對太陽射電進行觀測,得到的數(shù)據(jù)是強度隨輻射頻率和時間變化的圖譜,稱為射電動態(tài)頻譜圖。根據(jù)動態(tài)頻譜上的射電信號變化,我們可以找到頻率隨時間緩慢變化的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)基頻或二次諧頻會分裂成形態(tài)相似的兩個部分這一頻率分裂現(xiàn)象[5],該現(xiàn)象的特征參數(shù)對于理解II型射電暴的產(chǎn)生機制十分重要。近十幾年來,國內(nèi)外學(xué)者僅對少數(shù)精細結(jié)構(gòu)開展了特征提取研究,尚未涉及頻率分裂現(xiàn)象。所以,目前開展頻率分裂現(xiàn)象的特征自動提取研究將會對后續(xù)的太陽物理研究工作提供很大幫助。
中國科學(xué)院云南天文臺擁有兩臺具有了極高頻率分辨率和時間分辨率的全數(shù)字化太陽射電頻譜儀,能夠?qū)Ω鼮榫毢涂焖僮兓纳潆姳l(fā)及精細結(jié)構(gòu)進行觀測。本文實驗所采用的太陽射電頻譜數(shù)據(jù)來自云南天文臺撫仙湖太陽觀測基地的11米口徑米波太陽射電頻譜儀,這臺米波射電頻譜儀于2011年7月開始運行,工作頻段為70-700MHz,最高頻率分辨率為200kHz,優(yōu)于德國AIP的幾百KHz和澳大利亞Culgoora的1MHz分辨率[6],處于國際領(lǐng)先水平??紤]到時間分辨率和頻率分辨率之間的矛盾,該儀器日常工作的時間分辨率被設(shè)置為80ms,頻率分辨率設(shè)置為200kHz。工作時間為北京時間8:00-18:00。
表1 云南天文臺米波太陽射電頻譜儀的觀測參數(shù)
本文對原始頻譜中的通道做歸一化處理,消除強背景干擾和去除通道差異產(chǎn)生的噪聲[7],增強圖像中有用的信息。通道歸一化的原理是用頻譜圖像中每個通道的凈太陽射電爆發(fā)流量來除以該通道的平均值,以此減小各個通道之間的差異來實現(xiàn)噪聲消除。由于單通道的凈太陽射電爆發(fā)流量?fy(x)為原始通道的流量fy(x)減去該通道的平均值,所以單通道的歸一化公式為通過對y取全通道信息實現(xiàn)整幅圖像的通道歸一化。
大津法即最大類間方差法,能夠根據(jù)太陽射電頻譜圖像本身的灰度特性,將圖像分割為目標(biāo)前景與背景兩個部分。算法公式為分別為前景點數(shù)和背景點數(shù)所占圖像的比例,μa和μb分別為前景和背景的平均灰度。算法能自動選取最佳閾值T計算目標(biāo)前景與背景的最大類間總方差(T),實現(xiàn)二值分割。
經(jīng)典的精細結(jié)構(gòu)特征提取方法有連片搜索算法[9]、level set提取法[10]、小波變換分離法[11],這些方法僅適用于區(qū)域型精細結(jié)構(gòu)的提取。針對頻率分裂現(xiàn)象這類邊緣型的精細結(jié)構(gòu),本文提出了新的方法。
圖1 爆發(fā)區(qū)域確定實驗
高斯濾波是一種線性平滑濾波,常用于圖像的平滑處理。本文使用一個標(biāo)準(zhǔn)差為0.5,大小為5×5的卷積核對爆發(fā)區(qū)域進行二維高斯卷積平滑處理。對頻譜圖像進行平滑處理可以降低爆發(fā)區(qū)域噪聲,使爆發(fā)邊界更加清晰,增強精細結(jié)構(gòu)的特征。
為了消除不良干擾,保留爆發(fā)區(qū)域的重要結(jié)構(gòu)信息,需要對平滑處理后的爆發(fā)進行邊緣檢測。文獻[12]中,P.J.Zhang等人利用Canny邊緣檢測算法對來自NDA的射電數(shù)據(jù)進行邊緣檢測,得到很好的實驗效果。但是采用Canny算法對云南天文臺的數(shù)據(jù)進行處理會加強噪聲,而且檢測出的邊緣呈不閉合狀態(tài),其他傳統(tǒng)邊緣檢測方法的效果同樣不佳,因此本文根據(jù)邊緣指示函數(shù)[13]設(shè)計了適合的算法。
距離正則化水平集是以梯度流形式進行演化的,能夠有效處理分裂與合并這一類型的拓?fù)渥兓痆14],本文基于Li[14]等人提出的DRLSE算法進行編程實現(xiàn)并應(yīng)用在基于邊緣的主動輪廓模型上,實現(xiàn)距離正則化水平集方法自動提取精細結(jié)構(gòu)。針對頻率分裂現(xiàn)象的不規(guī)則性,該方法能很好的實現(xiàn)輪廓特征提取,而且演化過程非常穩(wěn)定。
(1)初始輪廓設(shè)置
太陽射電II型爆發(fā)發(fā)生在頻率較低的區(qū)域,起始頻率不會超過幾百兆赫茲。2007年White S針對觀測到的頻率分裂高頻事件報道了頻率分裂現(xiàn)象的基頻起始頻率為380MHz,二次諧頻的起始頻率為760MHz[6]這一研究成果,這為頻率分裂現(xiàn)象的觀測和進一步研究提供了一大依據(jù)。因此,結(jié)合本文實驗數(shù)據(jù)自身的頻段信息,將初始輪廓位置設(shè)置在380MHz到700MHz的范圍內(nèi)。
(2)頻率分裂現(xiàn)象的自動提取
距離正則化水平集方法依賴于圖像的梯度和能量函數(shù)實現(xiàn)演化功能,能量函數(shù)最小化運算驅(qū)動輪廓以梯度流形式進行演化,圖2(c)中的等高線即為能量值變化的體現(xiàn),初始輪廓會向著梯度較大的邊緣不斷進行能量最小化運算以接近邊緣。當(dāng)輪廓越來越接近邊緣時,算法會減緩輪廓的演化直至停止。本文根據(jù)經(jīng)驗將演化的迭代次數(shù)設(shè)為固定值1200,能夠滿足云南天文臺數(shù)據(jù)的精細結(jié)構(gòu)提取。
實 驗定義能量函數(shù)為E(?) =μRp(?) +Eext(?),其中μRp(?)為正則項,能夠保持演化的穩(wěn)定性。外部能量函數(shù)Eext(?)中包含邊緣項和區(qū)域向 :Eext(φ) =λLg(φ) +αA(φ),通過最小化外部能量函數(shù)能夠?qū)崿F(xiàn)演化功能。其中,區(qū)域項αA(φ)能夠確定演化方向,邊緣項λLg(φ)取最小值時說明輪廓已經(jīng)演化到目標(biāo)邊界處,演化停止。
(3)精細結(jié)構(gòu)提取結(jié)果及分析
實驗得到如圖2(d)中所示的提取結(jié)果,說明本文方法確實能有效實現(xiàn)頻率分裂現(xiàn)象的特征自動提取。由于本文數(shù)據(jù)中存在破壞了II型射電暴完整性的壞通道,導(dǎo)致輪廓在垂直方向上的演化最終停在射電暴的截斷處,但這并不影響最終的提取效果。假設(shè)沒有壞通道的干擾,輪廓在迭代1200次后同樣會在恰當(dāng)?shù)奈恢猛V?,最終提取出頻率分裂現(xiàn)象,而不會肆意演化。
圖2 頻率分裂現(xiàn)象提取實驗
精細結(jié)構(gòu)中的大部分特征參數(shù)只與其所在位置有關(guān)與爆發(fā)強度無關(guān)[15]。所以在提取出頻率分裂現(xiàn)象后,可以通過記錄頻率分裂現(xiàn)象的頻率最大值點和最小值點所在的通道和時間來計算精細結(jié)構(gòu)的持續(xù)時間和帶寬,然后利用線性回歸方法計算出頻率漂移率[16]。
本文的不足之處在于滿足研究要求的II型爆發(fā)數(shù)據(jù)不多,所以在提取出精細結(jié)構(gòu)之后只是簡單描述了可行的特征參數(shù)統(tǒng)計分析方法,希望以后在更多數(shù)據(jù)積累下,實現(xiàn)完善的特征參數(shù)統(tǒng)計分析工作。另外,II型射電爆發(fā)中還存在Herringbone精細結(jié)構(gòu),常常出現(xiàn)在頻率分裂現(xiàn)象的基頻和二次諧頻上,在本文的基礎(chǔ)上,可以開展相關(guān)研究。