基于改進混合高斯模型的太陽Hα圖像特征檢測

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

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

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

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

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

1 三幀差分法

幀間差分法是圖像處理過程中最常用的一種方法。該算法是將相鄰兩幀圖像之間做差來獲得運動區域。通過相鄰圖像的差值可以快速檢測出相鄰圖像中目標的運動范圍。該算法計算速度快，且對于光照變化不敏感。三幀差分法是取連續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)廣泛應用于數據分析和模式識別等領域，由Stauffer等人[18]在1999年提出。用多個高斯分布(正態分布)來表示圖像中t時刻每個像素點的灰度值Xi,t所構成的模型被稱作混合高斯模型，其函數表達式如式(4)所示，其中，P(Xt)表示每個像素點K個高斯分布的總概率；η(Xt,μi,t,∑i,t)表示像素點的概率密度；ωi,t表示第i個高斯分布t時刻的權重；μi,t表示第i個高斯分布t時刻的期望；∑i,t表示協方差矩陣；K表示高斯分布的個數，一般取值為3～5，本文取4。隨著K值的增大，模型所表示的場景就越復雜，但計算量也隨之增加。

(4)

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

(9)

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

3.1 圖像預處理

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

3.2 運動目標提取與檢測

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

3.3 太陽活動相關數據計算

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

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

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

表1 耀斑等級劃分標準

4 實驗結果與分析

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

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

圖1 算法比較圖

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

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

圖3 算法比較圖

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

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

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

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

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

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

圖6 GOES衛星給出的耀斑數據(Solarmonitor網站)

表2 3種算法對應的誤檢率

5 討論

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