沈發(fā)新，高冠男，汪 敏

(1. 中國科學(xué)院云南天文臺，云南 昆明 650216；2. 中國科學(xué)院大學(xué)，北京 100049；3. 中國科學(xué)院天體結(jié)構(gòu)與演化重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650216)

作為一個非?；钴S的天體，太陽上經(jīng)常會有劇烈活動或爆發(fā)現(xiàn)象(例如耀斑和日冕物質(zhì)拋射)。耀斑和日冕物質(zhì)拋射過程中產(chǎn)生的高能粒子和激波等在日地空間傳播，有可能到達(dá)地球，由此引發(fā)一系列日地物理效應(yīng)(例如太陽高能粒子事件、地磁暴等)，對人類的航天工程、通訊設(shè)施以及長距離輸油管道等都有非常重要的影響。因此，為了對空間天氣進(jìn)行有效的預(yù)報和預(yù)警，減少災(zāi)害性空間天氣對人類生產(chǎn)生活的破壞，深入細(xì)致地研究耀斑和日冕物質(zhì)拋射過程中產(chǎn)生的高能粒子和激波非常必要。

一般來說，天體射電信號十分微弱，但太陽射電暴的信號略強(qiáng)[1]。太陽射電爆發(fā)是指太陽活動期間在射電波段出現(xiàn)的劇烈且短促的流量增強(qiáng)現(xiàn)象，根據(jù)在射電動態(tài)頻譜上展現(xiàn)的不同形態(tài)，可以分為太陽射電I-V型射電暴。太陽射電Ⅲ型和Ⅱ型暴是目前研究比較多的兩種類型，分別是由耀斑和日冕物質(zhì)拋射期間產(chǎn)生的高能電子束流和激波引起的。

具體來說，耀斑和日冕物質(zhì)拋射期間的磁重聯(lián)加速電子形成高能電子，高能電子束流以亞相對論速度(約為0.2～0.6倍光速)沿磁力線快速運(yùn)動引起周圍等離子振蕩產(chǎn)生射電Ⅲ型暴，因此，射電Ⅲ型暴在射電動態(tài)頻譜上表現(xiàn)出快速的頻率漂移，是高能電子束流的最佳示蹤器。

日冕物質(zhì)拋射在日冕和行星際空間快速運(yùn)動時，當(dāng)日冕物質(zhì)拋射的速度超過本地的阿爾芬速度時，會產(chǎn)生日冕激波或行星際激波。通常認(rèn)為太陽射電Ⅱ型暴是由于激波在日冕中向外傳播引起等離子體振蕩產(chǎn)生朗繆波，再轉(zhuǎn)換為電磁波輻射并以本征等離子體頻率和二倍頻向外輻射。所以，射電Ⅱ型暴的觀測特征一般具有基頻和二次諧頻結(jié)構(gòu)，由于激波的運(yùn)動速度遠(yuǎn)低于高能電子束流，因此，太陽射電Ⅱ型暴在射電動態(tài)頻譜上通常表現(xiàn)為緩慢的頻率漂移。射電Ⅱ型暴可以視作日冕物質(zhì)拋射激波的最佳示蹤器。

目前，通常認(rèn)為太陽射電Ⅲ型和Ⅱ型暴的主要輻射機(jī)制為等離子體輻射，根據(jù)等離子體輻射的特點(diǎn)可以通過射電暴的信號頻率直接得到射電輻射源所在區(qū)域的密度，再根據(jù)日冕密度模型可以獲得輻射源所處的日冕高度，進(jìn)而獲得輻射源的運(yùn)動速度，例如日冕物質(zhì)拋射激波速度[1]。

現(xiàn)有的太陽射電觀測設(shè)備普遍有較高的頻率分辨率和時間分辨率，數(shù)據(jù)量較大，如果僅采用人工識別的方法，無法滿足實(shí)際需求。到目前為止, 世界上已經(jīng)發(fā)展了多種太陽射電暴識別算法。這些算法從不同角度處理觀測數(shù)據(jù)，如對流量數(shù)據(jù)使用閾值法確認(rèn)有無爆發(fā)[2]，對頻譜圖像數(shù)據(jù)使用霍夫變換和主動輪廓的方法確定爆發(fā)事件的頻率-時間關(guān)系[3]。

本文在這些識別算法的基礎(chǔ)上，使用概率霍夫變換的方法提取太陽射電Ⅲ型和Ⅱ型暴的起止時間、起止頻率、頻率漂移率等重要參數(shù)，提高了識別的準(zhǔn)確率和速度，并利用識別出的射電暴的物理參數(shù)結(jié)合相應(yīng)的物理模型估計了日冕物質(zhì)拋射激波的速度。這不但對空間天氣的預(yù)報和預(yù)警有重要的實(shí)用價值，而且對研究高能粒子束流和日冕物質(zhì)拋射在日冕和行星際空間傳播物理過程也具有重要的意義。

1 射電Ⅲ型暴的自動識別

射電Ⅲ型暴是一種最常見的太陽射電暴，是太陽劇烈活動后產(chǎn)生的非熱高能電子束在日冕和行星際空間運(yùn)動的最佳示蹤器[4]，同時也能最先反映太陽暴能量釋放和粒子加速等物理過程，是耀斑、日冕物質(zhì)拋射等太陽爆發(fā)事件的先兆。

由圖1可知，在不同的頻率范圍內(nèi)，太陽射電Ⅲ型暴的頻譜特征不同。例如，低于10 MHz，由于Ⅲ型射電暴輻射源所處的日冕密度隨著日冕高度的變化較慢，因此，10 MHz以下的射電Ⅲ型暴的頻率漂移速度遠(yuǎn)低于10 MHz以上的。而高于10 MHz的射電Ⅲ型暴具有類似直線的結(jié)構(gòu)，這使得本文使用概率霍夫變換算法進(jìn)行射電Ⅲ型暴識別成為可能。

圖1 Phoenix 4射電頻譜儀、南希10米波射電頻譜儀和WIND/WAVES射電探測器在2014年2月16日獲得的一例超寬帶射電Ⅲ型暴的綜合頻譜圖[5]

本文使用法國南希10米波天線陣(Nancay Decameter Array, NDA)的太陽射電觀測數(shù)據(jù)進(jìn)行Ⅲ型暴識別，其時間分辨率為1 s，觀測頻率范圍為10～80 MHz，其中均勻分布著400個通道。

為了識別觀測數(shù)據(jù)中包含的太陽射電暴，我們需要對觀測數(shù)據(jù)進(jìn)行處理。處理過程 (圖2) 分為3部分：(1)原始動態(tài)頻譜圖分段及背景去除；(2)動態(tài)頻譜圖前景提取及噪聲處理；(3)識別Ⅲ型暴事件并提取爆發(fā)信息。

圖2 Ⅲ型暴識別流程

1.1 太陽射電頻譜數(shù)據(jù)預(yù)處理

在太陽射電動態(tài)頻譜圖上，Ⅲ型暴表現(xiàn)為具有快速頻率漂移的信號，且不同Ⅲ型暴事件的相對強(qiáng)度有很大差異，而背景信號強(qiáng)度在幾小時內(nèi)逐漸變化。因此，我們把原始數(shù)據(jù)切分為若干小片段，這樣有助于去除背景效應(yīng)，同時也能盡量避免信號強(qiáng)的射電暴和弱的射電暴在同一個數(shù)據(jù)段中，導(dǎo)致弱的射電暴不能識別。經(jīng)過多次試驗(yàn)測試，我們得出每個片段為10 min的觀測數(shù)據(jù)可以有助于識別并提取射電暴的參數(shù)。

在南希10米波天線陣可觀測的頻率范圍內(nèi)，我們將每個數(shù)據(jù)段的時間長度設(shè)置為600 s，采樣間隔500 s，即兩個相鄰數(shù)據(jù)段的重疊時間為100 s。這種切分方式可以保證在某個片段中恰好切斷的爆發(fā)事件在下一個片段中可以識別。 因此，切分之后的數(shù)據(jù)段都是[600 × 250]的二維數(shù)組I(ti,fj)。

由于觀測設(shè)備不同通道之間存在明顯的增益差異，由此產(chǎn)生的通道效應(yīng)在頻譜圖上表現(xiàn)為高亮的橫條紋，這些橫條紋造成爆發(fā)結(jié)構(gòu)間斷，有時甚至遮蓋較弱的爆發(fā)事件，影響對太陽射電暴的識別，所以必須先消除觀測數(shù)據(jù)中的通道效應(yīng)。具體處理過程為在各通道的數(shù)據(jù)上減去該通道所有數(shù)據(jù)的均值，以去除背景，只保留來自太陽的射電信號，實(shí)現(xiàn)各通道數(shù)據(jù)歸一化，具體公式[3]為

(1)

其中，I(ti,fj)和I′(ti,fj)分別表示處理前和處理后的數(shù)據(jù)；a和b分別表示該數(shù)據(jù)段的開始時刻和結(jié)束時刻的索引。對數(shù)據(jù)段I中的每個值進(jìn)行計算，生成消除通道效應(yīng)后的I′。

1.2 Yen算法二值化及圖像腐蝕

對太陽射電動態(tài)頻譜圖進(jìn)行預(yù)處理后，我們還需要對圖像進(jìn)行二值化，將圖像灰度值置為0或255，經(jīng)過歸一化后將灰度值置為0或1，使目標(biāo)前景與背景區(qū)分開來，以消除背景的影響，有利于尋找目標(biāo)前景中的太陽射電Ⅲ型暴事件。

太陽射電動態(tài)頻譜圖二值化算法有很多種，但效果各有優(yōu)劣。本文使用多種自動確定閾值的算法對觀測數(shù)據(jù)進(jìn)行測試，生成二值化圖像，原始圖像與使用各種閾值法進(jìn)行二值化的圖像如圖3。

在圖3中，Original代表原始圖像，其余圖像代表各種算法(如Isodata, Li, Mean, Minimum, Otsu, Triangle和Yen)對原始圖像進(jìn)行二值化的輸出結(jié)果。其中最為典型的大津法(即Otsu法)[6]，也稱為最大類間方差法，是一種使用較為廣泛的自適應(yīng)閾值化方法。該算法通過不斷嘗試各種閾值，計算圖像的類間方差，選取類間方差最大時的閾值作為對該圖像進(jìn)行二值化的標(biāo)準(zhǔn)。此時，原圖像以該閾值為界限，將圖像中的所有像素分為背景像素和目標(biāo)像素兩大類。

圖3 各種自動確定閾值算法的二值化效果

從圖3的測試結(jié)果可以看出，只有Yen算法[7]能較好地確定閾值，在二值化后明確區(qū)分圖像前景與背景，并較好地保留前景信息，即爆發(fā)的信息，未與背景信息混淆。而其他算法的輸出結(jié)果未能準(zhǔn)確區(qū)分圖像前景與背景，導(dǎo)致爆發(fā)信息被干擾，難以識別。

Yen算法是一種自動多門限閾值法，能夠根據(jù)太陽射電頻譜圖像的灰度特性，將圖像分割為目標(biāo)前景和背景兩部分。與Otsu法類似，Yen算法能自動選取最佳閾值，實(shí)現(xiàn)圖像的二值分割，但其選擇最佳閾值的依據(jù)是最大熵準(zhǔn)則，基本思想是通過選擇閾值使背景和目標(biāo)共同提供的信息量最大。

如圖4，對太陽射電動態(tài)頻譜圖進(jìn)行二值化后，我們還需要對圖像進(jìn)行腐蝕(Erosion)操作，去除圖像中約每小時一次的定標(biāo)信號，以及某些在頻譜圖中占面積較小的瞬時干擾信號。

腐蝕是一種基本的形態(tài)學(xué)運(yùn)算，可以尋找圖像中的極小區(qū)域，將圖像中的高亮區(qū)域進(jìn)行縮減，在輸出的圖像中，高亮區(qū)域比原圖的更小。

從圖4可以看出，腐蝕算法較好地去除了頻譜圖中的定標(biāo)信號與瞬時干擾。同時，受腐蝕算法的影響，爆發(fā)事件結(jié)構(gòu)的邊緣也遭到腐蝕，但整體結(jié)構(gòu)與頻率漂移率不變，不影響后續(xù)識別。

圖4 經(jīng)過二值化與腐蝕后的圖像Fig.4 Image after binarization and erosion

1.3 使用概率霍夫變換算法識別太陽射電暴事件

霍夫變換是一種特征檢測方法，廣泛應(yīng)用于圖像分析、計算機(jī)視覺等領(lǐng)域。霍夫變換(圖5)可以在圖像中尋找某些圖形，如直線、圓和橢圓等。

圖5 霍夫變換原理[8]Fig.5 Hough transform principle[8]

常用的霍夫變換檢測直線的方法是利用

ρ=xcosα+ysinα

(2)

在圖像空間和參數(shù)空間之間,將圖像中的曲線轉(zhuǎn)換為參數(shù)空間的一點(diǎn),建立對偶變換。(2)式中，ρ為極徑；α為極角，取0～180°；x為像素點(diǎn)相對圖像原點(diǎn)的行坐標(biāo)；y為像素點(diǎn)相對圖像原點(diǎn)的列坐標(biāo)。

如圖6，輸入圖像從圖像空間轉(zhuǎn)換到參數(shù)空間后，在參數(shù)空間存在3個明顯的亮點(diǎn)，而參數(shù)空間的亮點(diǎn)代表輸入圖像中的直線結(jié)構(gòu)。利用霍夫變換的直線檢測功能，我們可以在太陽射電頻譜圖像中尋找擁有類似直線結(jié)構(gòu)的太陽射電暴。

圖6 霍夫變換的效果Fig.6 The effect of Hough transform

本文使用概率霍夫變換算法解析曲線的幾何特性，盡量避免對圖像中的每個非零像素進(jìn)行計算量大的投票過程，而是從點(diǎn)集中隨機(jī)選取一個像素點(diǎn)進(jìn)行分析。

概率霍夫變換算法具體步驟為(1)每個區(qū)間對應(yīng)一個累加器A(r,θ)，從每個區(qū)間的點(diǎn)集中隨機(jī)選取一個像素點(diǎn)，在參數(shù)空間計算θ對應(yīng)的r值，并在對應(yīng)的累計器A(r,θ)上加1；(2)從點(diǎn)集中刪除該點(diǎn)并更新累加器；(3)若累加器的值大于設(shè)定的閾值，則表示已檢測到直線，然后刪除該直線上所有的點(diǎn)；(4)重復(fù)以上步驟，直至點(diǎn)集為空。

根據(jù)太陽射電Ⅲ型暴的頻譜信號特征，爆發(fā)事件在太陽射電動態(tài)頻譜圖中近似一條傾斜的直線。直線的傾斜程度與太陽射電Ⅲ型暴的頻率漂移率相關(guān)。

我們對太陽射電動態(tài)頻譜圖使用概率霍夫變換尋找Ⅲ型暴。如果在動態(tài)頻譜數(shù)據(jù)段中存在符合約束條件的直線，則認(rèn)為該動態(tài)頻譜數(shù)據(jù)段中可能存在Ⅲ型暴事件。但是，在數(shù)據(jù)段中也可能存在其他類型的射電暴，例如太陽射電Ⅱ型暴在頻譜上也可能具有類似直線的圖像特征。所以需要根據(jù)識別的頻率漂移率來限制和過濾其他類型的爆發(fā)事件。

動態(tài)頻譜數(shù)據(jù)段經(jīng)過二值化、腐蝕處理后，對每一幀二值化后的圖像進(jìn)行概率霍夫變換，檢測圖像中的直線結(jié)構(gòu)。同時，已知在南希10米波天線陣觀測數(shù)據(jù)的頻率范圍內(nèi)，太陽射電Ⅲ型暴的頻率漂移率約為1～10 MHz/s，所以我們需要對檢出的直線進(jìn)行篩選，只留下頻率漂移率符合Ⅲ型暴的直線。若經(jīng)過篩選后，數(shù)據(jù)段中仍存在一條或多條直線，則確定該數(shù)據(jù)段中至少包含一個事件。

如圖7，圖中有直線特征的結(jié)構(gòu)被識別并標(biāo)紅，此處很可能是一個Ⅲ型暴事件，同時還可以獲取線段兩端的坐標(biāo)值，即爆發(fā)事件的起始時間、起始頻率、終止時間和終止頻率。圖8中事件的起始頻率為50 MHz，終止頻率為27 MHz，持續(xù)時間為15 s，頻率漂移速率為-1.5 MHz/s，負(fù)值表示Ⅲ型暴從高頻向低頻漂移。

1.4 南希10米波天線陣數(shù)據(jù)中Ⅲ型暴的識別結(jié)果

以上識別算法可以高效地在大于10 MHz的頻譜數(shù)據(jù)中識別形態(tài)類似直線的Ⅲ型暴事件，并獲得起止頻率、起止時間和頻率漂移率等信息，為研究爆發(fā)事件提供支持。

如圖8，圖8(a)中爆發(fā)事件的起始頻率為67 MHz，終止頻率為23 MHz，持續(xù)時間為17 s，頻率漂移率為-2.6 MHz/s。圖8(b)中爆發(fā)事件的起始頻率為37.3 MHz，終止頻率為28 MHz，持續(xù)時間為6 s，頻率漂移率為1.55 MHz/s。以上兩例爆發(fā)事件的識別結(jié)果證明，本文的識別算法可以識別從各頻率開始的爆發(fā)事件。圖8(c)中爆發(fā)事件的起始頻率為58 MHz，終止頻率為46 MHz，持續(xù)時間為6 s，頻率漂移率為-2.1 MHz/s。圖8(d)中爆發(fā)事件的起始頻率為48 MHz，終止頻率為26 MHz，持續(xù)時間為18 s，頻率漂移率為-1.22 MHz/s。以上兩例爆發(fā)事件的識別結(jié)果證明，本文的識別算法可以識別各種長度的爆發(fā)事件。

圖8 識別出的不同參數(shù)的爆發(fā)事件Fig.8 Burst events with different parameters identified

但由于南希10米波天線陣觀測的頻率范圍有限，且在觀測范圍內(nèi)的低頻頻段中存在較強(qiáng)持續(xù)不斷的干擾，如果爆發(fā)事件的終止頻率較低，這些干擾信號可能影響算法對爆發(fā)事件終止頻率的識別。如圖9，算法識別的爆發(fā)事件的終止頻率與干擾信號所在頻率緊密結(jié)合，這樣在識別爆發(fā)事件終止頻率時有一定的誤差。

圖9 識別出的低頻爆發(fā)事件Fig.9 Low-frequency burst event identified

2 射電Ⅱ型暴的自動識別

從太陽射電觀測數(shù)據(jù)中識別太陽射電Ⅱ型暴的方法與識別Ⅲ型暴的方法類似，需要對數(shù)據(jù)進(jìn)行如圖2的處理過程。但在使用概率霍夫變換算法識別頻譜圖中的直線結(jié)構(gòu)及起止點(diǎn)坐標(biāo)后，需要根據(jù)Ⅱ型暴的物理特征進(jìn)行篩選。例如Ⅱ型暴的基波所在的起始頻率通常不高(<300 MHz)，頻率漂移速率較為緩慢(<1 MHz/s)。

Ne=N0104.32Rs/R，

(3)

其中，N0=4.2×104cm-3；Rs為太陽半徑。

基于以上自動提取算法和日冕密度模型，我們對于2012年11月13日云南天文臺米波太陽射電頻譜儀的觀測數(shù)據(jù)進(jìn)行處理，識別其中的Ⅱ型暴，提取起止頻率、起止時間、頻率漂移率等參數(shù)，并計算所對應(yīng)的日冕物質(zhì)拋射的特征參數(shù)(如圖10)。

圖10 利用觀測數(shù)據(jù)提取日冕物質(zhì)拋射特征參數(shù)Fig.10 Extracting CME feature parameters from observation data

通過概率霍夫變換算法得到該Ⅱ型射電暴的物理參數(shù)，并得出對應(yīng)的日冕物質(zhì)拋射在02:02:20 UT-02:05:44 UT的速度為870 km/s，而后，LASCO C2日冕儀在02:24:06 UT發(fā)現(xiàn)該日冕物質(zhì)拋射，并獲得該日冕物質(zhì)拋射的速度為851 km/s(數(shù)據(jù)來源：https://cdaw.gsfc.nasa.gov/CME_list/UNIVERSAL/2012_11/univ2012_11.html)。通過對比發(fā)現(xiàn)，這兩種方法給出的日冕物質(zhì)拋射速度非常接近，因此，我們認(rèn)為利用概率霍夫變換算法自動識別的參數(shù)是合理可信的。

3 總結(jié)與展望

本文提出了一種基于概率霍夫變換的太陽射電Ⅲ型和Ⅱ型暴識別算法。首先對太陽頻譜圖像進(jìn)行預(yù)處理，多次試驗(yàn)測試選擇合適的時間長度對原始觀測數(shù)據(jù)進(jìn)行合理切分；對各通道進(jìn)行均值化，消除通道不均勻的干擾；對圖像進(jìn)行二值化，提取前景信息；對圖像進(jìn)行腐蝕，消除一些窄帶射頻干擾信號。然后使用概率霍夫變換算法尋找頻譜圖像中的直線結(jié)構(gòu)，獲?、笮秃廷蛐捅┦录钠鹬箷r間、起止頻率和頻率漂移率等信息，并根據(jù)這些信息進(jìn)一步篩選頻率漂移率符合Ⅲ型和Ⅱ型暴特征的直線結(jié)構(gòu)。

本文還利用該方法從2012年11月13日云南天文臺的觀測數(shù)據(jù)中自動提取了Ⅱ型暴的物理參數(shù)，并結(jié)合日冕密度模型獲得日冕物質(zhì)拋射的特征參數(shù)。這些參數(shù)對空間天氣的預(yù)報預(yù)警具有實(shí)用價值。該方法雖然可以識別典型的Ⅱ型暴，但由于Ⅱ型暴容易與其他類型射電暴同時發(fā)生，信號有混疊，未來我們將在更多的Ⅱ型暴數(shù)據(jù)上對該算法進(jìn)行檢驗(yàn)和改進(jìn)，以達(dá)到更好的、對多種數(shù)據(jù)來源普遍適用的識別效果。

本文結(jié)果表明，該方法能有效地識別太陽射電Ⅲ型暴和Ⅱ型暴，獲取射電暴的起止頻率、起止時間和頻率漂移率等信息，并能快速地對長期、大量的觀測數(shù)據(jù)進(jìn)行識別和分析。同時，概率霍夫變換算法比標(biāo)準(zhǔn)霍夫變換算法的計算量更小，計算過程更短，更適合在觀測過程中實(shí)時識別太陽射電Ⅲ型和Ⅱ型暴，提取觀測參數(shù)得出日冕物質(zhì)拋射激波速度等關(guān)鍵物理信息，實(shí)現(xiàn)空間天氣預(yù)報預(yù)警，維護(hù)相關(guān)人員、設(shè)備的安全。