梁宇峰，白先勇,2,3，馮 松

(1. 昆明理工大學(xué) 云南省計算機技術(shù)應(yīng)用重點實驗室/信息工程與自動化學(xué)院,云南 昆明 650500；2.中國科學(xué)院 太陽活動重點實驗室,北京100012；3.中國科學(xué)院 國家天文臺,北京 100049)

太陽活動是太陽大氣中局部區(qū)域的各種不同活動現(xiàn)象的總稱。當(dāng)太陽活動增強時會突然釋放出巨大的能量，同時拋射出不同能量的粒子對地球空間環(huán)境造成很大的影響。隨著計算機技術(shù)的快速發(fā)展，利用運動目標檢測技術(shù)對太陽圖像進行檢測，從而有效地監(jiān)測太陽活動是目前天文圖像處理的研究熱點之一。由于受到大氣吸收的影響，目前在地面可以觀測到太陽的光球?qū)雍蜕驅(qū)?，而日冕層的觀測還僅限于日面邊緣。Hα圖像可以充分顯示色球豐富結(jié)構(gòu)和各類色球活動現(xiàn)象，是最有效的色球觀測譜線，國內(nèi)的觀測設(shè)備如中國科學(xué)院云南天文臺全日面Hα色球望遠鏡[1]、懷柔太陽觀測基地的全日面Hα望遠鏡[2-3]以及光學(xué)和紅外太陽爆發(fā)監(jiān)測望遠鏡[4]。國外的觀測設(shè)備如美國國立太陽天文臺的太陽全球振蕩監(jiān)測網(wǎng)(Global Oscillation Network Group，GONG)[5-6]、全球高分辨率Hα網(wǎng)[7](Global High-Resolution H-alpha Network)都使用該譜線觀測色球活動。檢測太陽Hα圖像的特征也就是檢測太陽耀斑、日珥、黑子和暗條等太陽活動。

近年來，學(xué)者們對Hα圖像進行實驗并開發(fā)了一系列的檢測算法和程序。針對太陽暗條的檢測方法主要可以分為兩類。(1)是基于圖像形態(tài)處理方法。Hao等[8]采用形態(tài)學(xué)操作的方法，針對Mauna Loa Solar天文臺觀測到的Hα全日面太陽圖像進行實驗，識別暗條的結(jié)果較好。Fuller[9]采用了經(jīng)典的區(qū)域生長方法對法國Meudon天文臺的Hα圖像進行驗證，并得到了較好的結(jié)果，但是圖像局部差異容易導(dǎo)致算法的不穩(wěn)定性。Gao等[10]采用全局閾值和區(qū)域生長方法對美國Big Bear Solar天文臺的Hα圖像進行驗證，此法簡單高效，但檢測不到邊緣暗條，且易漏掉形狀較小暗條，黑子也極有可能混雜其中。(2)基于統(tǒng)計機器學(xué)習(xí)的方法。Zharkova[11-12]利用神經(jīng)網(wǎng)絡(luò)方法來識別暗條，并獲得了較好的效果。Qu[13]在文獻中采用支持向量機的方法來區(qū)分暗條和黑子。太陽耀斑的檢測也取得了較快的發(fā)展，Qu等[14]在2003年利用支持向量機分類器的方法對耀斑進行識別，取得了良好的效果。Borda等[15]利用神經(jīng)網(wǎng)絡(luò)算法實現(xiàn)了耀斑的自動檢測。P?tzi等[16]針對Kanzelh?he天文臺的Hα圖像做臨邊昏暗改正，并把存在的太陽活動做成高斯直方圖，這樣就有個灰度從而得到耀斑。從上述可以看出，如何分析海量的觀測數(shù)據(jù)，并開發(fā)更有效地算法來檢測Hα圖像是當(dāng)今的研究熱點。

本文將改進的混合高斯模型方法引入到太陽活動的檢測中。預(yù)處理圖像后，利用三幀差分法先得到的一幅背景圖像，再將該圖像更新到混合高斯模型的背景模型中，對正在爆發(fā)的太陽活動進行檢測并識別。通過對觀測圖像的處理表明：本文算法同時檢測出耀斑和暗條爆發(fā)等太陽活動，并計算出耀斑爆發(fā)時間和級別。通過對中國科學(xué)院國家天文臺懷柔太陽觀測基地和GONG觀測的3組數(shù)據(jù)進行實驗，結(jié)果證明本算法和現(xiàn)有的算法相比能有效的檢測出太陽耀斑、暗條等活動并計算出耀斑級別和時間。改善了混合高斯模型固有的光照變化檢測不佳的問題，提高了檢測耀斑、暗條等太陽活動的準確度。

1 三幀差分法

幀間差分法是圖像處理過程中最常用的一種方法。該算法是將相鄰兩幀圖像之間做差來獲得運動區(qū)域。通過相鄰圖像的差值可以快速檢測出相鄰圖像中目標的運動范圍。該算法計算速度快，且對于光照變化不敏感。三幀差分法是取連續(xù)3幀圖像Ii-1(x,y)、Ii(x,y)和Ii+1(x,y)。分別計算兩幀間的差值，如式(1)和式(2)所示[17]，選取合適的閾值T得到b(i,i-1)和b(i+1,i)。再對b(i,i-1)和b(i+1,i)做邏輯“與”運算得到運動目標Bi，如式(3)所示。

(1)

(2)

(3)

2 混合高斯模型方法

混合高斯模型(Gaussian Mixture Models，GMMs)廣泛應(yīng)用于數(shù)據(jù)分析和模式識別等領(lǐng)域，由Stauffer等人[18]在1999年提出。用多個高斯分布(正態(tài)分布)來表示圖像中t時刻每個像素點的灰度值Xi,t所構(gòu)成的模型被稱作混合高斯模型，其函數(shù)表達式如式(4)所示，其中，P(Xt)表示每個像素點K個高斯分布的總概率；η(Xt,μi,t,∑i,t)表示像素點的概率密度；ωi,t表示第i個高斯分布t時刻的權(quán)重；μi,t表示第i個高斯分布t時刻的期望；∑i,t表示協(xié)方差矩陣；K表示高斯分布的個數(shù)，一般取值為3～5，本文取4。隨著K值的增大，模型所表示的場景就越復(fù)雜，但計算量也隨之增加。

(4)

若Xi,t不匹配，當(dāng)k

ωi,t=(1-α)ωi,t-1+α(Mi,t)

(5)

μi,t=(1-ρ)μi,t-1+ρXi,t

(6)

(7)

ρ=αη(Xt|μi,t,σi,t)

(8)

(2)混合高斯模型的前景檢測。首先將每個像素的K個高斯分布根據(jù)ωi,t/σi,t比值大小進行排序。再從K個高斯分布中選取B個高斯分布作為背景模型。其中式(9)給出了B的取值。T表示背景閾值。將像素值Xi,t與B個高斯分布進行匹配，若存在匹配，則該像素點為背景點，否則該像素點被檢測為前景點。對圖像序列中的每個像素點重復(fù)執(zhí)行以上步驟，完成圖像的前景分割過程得到Hα圖像中的運動物質(zhì)。

(9)

3 太陽Hα圖像特征檢測方法過程

3.1 圖像預(yù)處理

使用改進混合高斯模型檢測太陽活動之前，先對Hα圖像進行預(yù)處理。預(yù)處理包括圖像的標準化和圖像復(fù)原兩個步驟。本文首先對Hα圖像進行標準差歸一化。使圖像符合標準高斯分布(即均值為0，標準差為1)，保證圖像的強度統(tǒng)一。在Hα圖像獲取時，受到天氣條件等各種隨機噪聲的影響，導(dǎo)致圖像存在信噪比低、圖像運動模糊等問題。因此要提高圖像質(zhì)量，必須對其進行濾波處理。維納濾波的原理是假設(shè)圖像信號可以近似看成平穩(wěn)隨機過程的前提下，使輸入圖像和復(fù)原圖像之間的均方誤差達到最小的準則函數(shù)來實現(xiàn)圖像復(fù)原。本實驗數(shù)據(jù)的模糊主要是由于天氣原因(云，陰天)和望遠鏡和場景之間的均勻線性運動產(chǎn)生的，因此再利用MATLAB工具箱[21]中的fspecial函數(shù)建模并用維納濾波方法(deconvwnr函數(shù))來對圖像進行平滑處理，恢復(fù)運動模糊。

3.2 運動目標提取與檢測

預(yù)處理結(jié)束后，采用改進混合高斯模型進行圖像的目標(太陽活動)提取。首先確定高斯分布的個數(shù)K。本文分別統(tǒng)計了3組數(shù)據(jù)中多個位置的背景像素的變化，發(fā)現(xiàn)3組數(shù)據(jù)集中大部分的模擬背景像素分別符合4個高斯分布的疊加，因此本文采用4個高斯分布作為數(shù)據(jù)集的模擬高斯的個數(shù)。接著初始化混合高斯背景模型，經(jīng)過實驗，本文采用圖像序列第1幅Hα圖像做初始背景，以第1幅圖像同位置的像素值作為對應(yīng)的均值，同時方差取比較大的值，第1個高斯分布取較大的權(quán)重(本文取0.994)，其他都取較小值來進行初始化。從理論上來說，混合高斯背景建模應(yīng)該是一種較為完美的背景分割方法。但在發(fā)生光線突變時，導(dǎo)致其檢測效果下降，為了解決該問題本文將三幀差分法對于光照變化不敏感的優(yōu)勢和混合高斯模型不容易出現(xiàn)“空洞”現(xiàn)象的優(yōu)勢相結(jié)合，克服了光照漸變帶來的影響。初始化后讀入數(shù)據(jù)，利用三幀差分法先得到運動區(qū)域F1，用當(dāng)前幀與F1做差得到F2，將F2作為背景區(qū)域更新到混合高斯模型的主背景分布中建立圖像序列中每個時刻圖像的動態(tài)背景模型[22]。符合高斯分布后進行參數(shù)更新。用觀測圖像減去背景得到前景圖像，得到的前景圖像就是正在爆發(fā)的太陽活動。雖然改善了光照突變帶來的影響，但在提取出來的前景圖像中依然會存在面積較小的干擾點。這些點會對后續(xù)的檢測造成一定的影響。由于暗條、耀斑面積通常大于3×108km2[24]，所以小于該值的面積一般為小耀斑，引起顯著日地效應(yīng)的概率不大，所以不是本文的關(guān)注對象。(將3×108km2換算到實驗數(shù)據(jù)中圖像上像素的大小，經(jīng)過計算在數(shù)據(jù)集Ⅰ上面積為376，數(shù)據(jù)集Ⅱ上為468)。因此本文算法利用MATLAB工具[21]箱中的bwareaopen函數(shù)清除掉這些區(qū)域，最后完成檢測操作。

3.3 太陽活動相關(guān)數(shù)據(jù)計算

(1)判斷太陽活動類型。由于暗條的平均亮度要低于太陽日面的平均亮度，而耀斑的平均亮度要高于太陽日面的平均亮度。因此檢測出太陽活動后標記在原圖像上。計算該區(qū)域強度的平均值和整個圖像強度的平均值，如果該區(qū)域強度平均值小于整個圖像的強度平均值，則該區(qū)域為暗條活動，反之為耀斑活動；

(2)計算耀斑活動時間。耀斑的開始時刻即是程序識別出耀斑的時刻；峰值時刻即是耀斑面積最大時刻；耀斑面積回到耀斑前狀態(tài)即為結(jié)束時刻。因此利用MATLAB軟件工具[21]箱中的regionprops函數(shù)分別計算每一幀耀斑(前景圖像)的面積來確定耀斑的開始、峰值和結(jié)束時刻；

(3)計算耀斑級別。傳統(tǒng)上，根據(jù)Hα耀斑極大時經(jīng)投影改正后的面積大小，把耀斑分為4個等級[24]。耀斑等級劃分標準[24]如表1所示(R⊙代表太陽半徑)。因此將識別的耀斑面積投影改正到日面中心的面積，再換算成平方度對照表1得到耀斑級別。

表1 耀斑等級劃分標準

4 實驗結(jié)果與分析

在Hα的色球觀測圖像中的太陽活動特征為耀斑、暗條和日珥。耀斑在日面和日面邊緣均能發(fā)生，暗條是發(fā)生在日面上的，在日面邊緣的叫日珥，因此挑選出的實驗數(shù)據(jù)集均包含上述特征。本研究選用了來做兩個不同觀測設(shè)備在Hα波段上的觀測資料。數(shù)據(jù)集 I 和 Ⅱ 來自中國科學(xué)院國家天文臺懷柔太陽觀測基地全日面光學(xué)和磁場監(jiān)測系統(tǒng)觀測的Hα色球資料。其中懷柔觀測的數(shù)據(jù)所采用的望遠鏡的主要參數(shù)：口徑20 cm，有效焦距180 cm，配備半寬0.25 ?的Lyot雙折射濾光器，濾光器放在準直光路中，探測器有效像元數(shù)為2 712×2 712[23-24]，像元分辨率為0.863 39″。數(shù)據(jù)集 Ⅲ 是來自Gong的Hα觀測資料數(shù)據(jù)集，其有效像元數(shù)為2 048×2 048。數(shù)據(jù)集 I 包含日面和日面邊緣的耀斑和暗條爆發(fā)，也包含1個日面上的小耀斑。數(shù)據(jù)集 Ⅱ 包含日面上的大耀斑。數(shù)據(jù)集 Ⅲ 包含日面上的中等強度耀斑和暗條爆發(fā)。所有數(shù)據(jù)集的時間分辨率均為60 s。數(shù)據(jù)集 I 的觀測時間是2013/05/12/ 23：01 UT到2013/05/13/ 04：22 UT共255幀圖像。數(shù)據(jù)集 II 是2013/05/14/ 22：42 UT到2013/05/15/ 03：05 UT共135幀圖像。數(shù)據(jù)集 Ⅲ 是從2015/11/04/ 03：10 UT到2015/11/04/ 08：48 UT共309幀圖像。

首先對數(shù)據(jù)集 I 進行實驗。預(yù)處理結(jié)束后設(shè)置三幀差分法的閾值T=27。設(shè)置GMMs的初始參數(shù)：采用4個高斯分布來模擬背景，權(quán)值第1位參數(shù)取0.994, 其余取等權(quán)0.002。方差為25，學(xué)習(xí)率α為0.005。經(jīng)過多次實驗調(diào)整，設(shè)置背景分布閾值T為0.7，均值第1位參數(shù)取第1幀圖像同位置的像素灰度值，其余均為0。經(jīng)過實驗得到結(jié)果。圖1第1行從左至右分別表示的是數(shù)據(jù)集 I 中的第146幀(2013/05/13/ 02：26UT)和147幀(2013/05/13/ 02：27UT)的部分原始Hα圖像，正在爆發(fā)的耀斑活動由圓圈出；圖1第2行表示三幀差分法對以上兩幀檢測出的前景圖像；圖1第3行表示傳統(tǒng)GMMs對以上兩幀檢測出的前景圖像；圖1第4行表示本文算法對以上兩幀檢測出的前景圖像。由圖1第2行可知，三幀差分得到的結(jié)果會有“空洞”現(xiàn)象；由圖1第3行可知，雖然檢測出了耀斑活動，但會有很多噪聲點對計算耀斑數(shù)據(jù)造成影響，導(dǎo)致識別率降低；由圖1第4行可知，本文算法可以較好地檢出耀斑區(qū)域。

圖1 算法比較圖

根據(jù)檢測結(jié)果計算得到結(jié)論：在日面邊緣位置(819″，517″)處編號為NOAA AR 11748的活動區(qū)有太陽耀斑爆發(fā)。根據(jù)耀斑面積隨時間的變化計算出耀斑的開始時刻為2013/05/13/ 02：09 UT、峰值時刻為2013/05/13/ 02：32 UT，結(jié)束時刻為2013/05/13/ 03：16 UT。經(jīng)過計算得到數(shù)據(jù)集Ⅰ的平方度約為9.673，查詢表1可知該耀斑級別為2。最后將爆發(fā)的耀斑活動位置標記在原始Hα圖像上。通過查詢Solarmonitor網(wǎng)站，可以得到地球靜止環(huán)境業(yè)務(wù)衛(wèi)星(Geostationary Operational Environmental Satellites，GOES)給出的耀斑數(shù)據(jù)，如圖6第1行。對比兩組時間可知，本文的檢測結(jié)果和GOES衛(wèi)星給出的時間段大致相同(時間上的偏差由于地面觀測時間分辨率設(shè)置的不同)。圖2中紅色部分表示的是耀斑活動(峰值時刻)在Hα圖像中的位置。

圖2 耀斑活動在原Hα圖像中的位置

圖3 算法比較圖

再對數(shù)據(jù)集 Ⅱ 進行實驗，三幀差分法的閾值T=30。設(shè)置GMMs的初始參數(shù)：采用4個高斯分布來模擬背景，權(quán)值第1位參數(shù)取0.994, 其余取等權(quán)0.002。方差為20，學(xué)習(xí)率α為0.006。背景分布閾值T為0.65，均值第1位參數(shù)取第1幀圖像同位置的像素灰度值，其余均為0。經(jīng)過實驗得到結(jié)果。圖3第1行左邊表示數(shù)據(jù)集 Ⅱ 中的第73幀(2013/05/15/ 01：51UT)的部分原始Hα圖像，正在爆發(fā)的耀斑活動由圓圈出；圖3第1行右邊表示三幀差分法對該幀檢測出的前景圖像；圖3第2行左邊表示傳統(tǒng)GMMs對該幀檢測出的前景圖像；圖3第2行右邊表示本文算法對該幀檢測出的前景圖像。由圖3可知，三幀差分法和傳統(tǒng)GMMs會帶來識別率降低的問題。本文算法可以很好的避免光照突變帶來的誤檢測現(xiàn)象，并檢測出耀斑區(qū)域。

根據(jù)檢測結(jié)果計算得到結(jié)論：在日面邊緣位置(739″，443″)處編號為NOAA AR 11748的活動區(qū)有太陽耀斑爆發(fā)。根據(jù)耀斑面積隨時間的變化計算出耀斑的開始時刻為2013/05/15/ 01：36 UT、峰值時刻為2013/05/15/ 01：54 UT，結(jié)束時刻為2013/05/15/ 02：16 UT。數(shù)據(jù)集 Ⅱ 的平方度約為6.521，查詢表1可知該耀斑級別為2。查詢Solarmonitor網(wǎng)站，得到GOES衛(wèi)星給出的耀斑數(shù)據(jù)，如圖6第2行。對比兩組時間可知，本文的檢測結(jié)果和GOES衛(wèi)星給出的時間段大致相同(時間上的偏差由于地面觀測時間分辨率設(shè)置的不同)。圖4中亮點部分表示的是耀斑活動(峰值時刻)在Hα圖像中的位置。

為驗證該算法的有效性，最后對數(shù)據(jù)集 Ⅲ 進行相同操作，并得到實驗結(jié)果。在2015/11/04/ 03：46 UT，日面位置為(167″，34″)處編號為NOAA AR 12443的活動區(qū)有劇烈的暗條活動發(fā)生。同時伴隨著1個較小級別的太陽耀斑爆發(fā)，經(jīng)過計算得到該耀斑的級別為S級。在日面位置(921″，256″)處編號為NOAA AR 12445的活動區(qū)有太陽耀斑爆發(fā)。耀斑的開始時刻為2015/11/04/ 03：20 UT、峰值時刻為2015/11/04/ 03：24 UT，結(jié)束時刻為2015/11/04/ 04：27 UT，耀斑級別為1。圖5中亮點部分表示的是耀斑活動(峰值時刻)在Hα圖像中的位置，灰條部分表示暗條活動和伴隨它發(fā)生的小級別耀斑活動在Hα圖像中的位置。GOES衛(wèi)星給出的耀斑數(shù)據(jù)如圖6第3行。本文的檢測結(jié)果和GOES衛(wèi)星給出的時間段大致相同。

圖4 耀斑活動在原Hα圖像中的位置

根據(jù)3組數(shù)據(jù)和GOES衛(wèi)星的數(shù)據(jù)對比可以驗證本文算法的有效行。表2是3種方法對耀斑檢測后計算得到的誤檢率。根據(jù)圖1、圖3和表2組可知，本文成功檢測耀斑和暗條的基礎(chǔ)上，改善了混合高斯模型固有的光照變化檢測不佳的問題，提高了檢測耀斑等太陽活動的準確度。

圖5 耀斑和暗條活動在原Hα圖像中的位置

圖6 GOES衛(wèi)星給出的耀斑數(shù)據(jù)(Solarmonitor網(wǎng)站)

表2 3種算法對應(yīng)的誤檢率

5 討論

太陽活動和地球空間環(huán)境有著密不可分的關(guān)系，對其活動爆發(fā)進行識別和跟蹤有著重要的意義。本文將改進的混合高斯模型引入到太陽活動檢測中，提出一種利用混合高斯背景差分模型建立動態(tài)背景來檢測Hα圖像上正在爆發(fā)的太陽活動識別算法。本文分別對3組來自中國科學(xué)院國家天文臺懷柔太陽觀測基地和GONG的數(shù)據(jù)進行實驗。首先對圖像進行歸一化和維納濾波處理，再利用三幀差法快速得到背景圖像，更新到混合高斯模型主背景分布中建立穩(wěn)定的Hα圖像動態(tài)背景，實現(xiàn)運動目標(太陽活動)的提取。本文算法在成功識別太陽耀斑、暗條活動，并計算出了耀斑、暗條活動的日面位置(邊緣、日面)。再根據(jù)耀斑面積的變化計算出耀斑的級別(C、M、X級)的同時，還提高了檢測到準確率。本算法的耀斑檢測結(jié)果和GOES衛(wèi)星給出的耀斑級別和時間一致，從而驗證了算法的有效性。相比GOES衛(wèi)星對整個日面的流量監(jiān)測，本算法不僅可以從Hα序列圖像的識別結(jié)果中給出耀斑級別，開始、峰值和結(jié)束時刻，還可以給出耀斑、暗條等太陽活動特征對應(yīng)的日面位置，對基于我國現(xiàn)有地基觀測的Hα色球數(shù)據(jù)開展太陽活動實時監(jiān)測預(yù)警具有重要意義。此外，GOES主要是反應(yīng)太陽日冕活動現(xiàn)象，相比較而言，Hα圖像反應(yīng)的則是更低層次即色球的物理過程，能更好地反應(yīng)太陽活動的源區(qū)特性，因此對于太陽物理有著重要意義。雖然改進混合高斯模型算法的會導(dǎo)致計算量變大，運算時間變長，但是卻提高了探測耀斑的準確度。對于耀斑而言提高其檢測的準確度比檢測速度對監(jiān)測預(yù)警來說有著更重要的作用。