基于支持向量機(jī)的太陽磁場活動區(qū)極性反轉(zhuǎn)線位置的探測?

王 強(qiáng) 鄭 勝 鄧元勇 黃 宇黎 輝 甘為群 蘇江濤,3

(1三峽大學(xué)計算機(jī)與信息學(xué)院宜昌443002)

(2中國科學(xué)院國家天文臺北京100101)

(3中國科學(xué)院大學(xué)北京100049)

(4中國科學(xué)院紫金山天文臺南京210023)

1 引言

太陽的活動現(xiàn)象十分豐富,其中太陽耀斑和日冕物質(zhì)拋射是兩類最劇烈的爆發(fā)現(xiàn)象.太陽耀斑是局部區(qū)域突然快速釋放大量能量的過程.耀斑爆發(fā)時釋放的大量高能離子對地球空間環(huán)境破壞性較大,尤其是影響航空航天與導(dǎo)航通信行業(yè).對太陽活動現(xiàn)象的研究是天文研究領(lǐng)域的一個重要方向.太陽活動區(qū)磁場極性反轉(zhuǎn)線的位置,長期以來一直是太陽活動觀測和理論研究的核心特征.

為了理解太陽耀斑的物理機(jī)制和開展研究耀斑預(yù)測的方法,最重要的是找到與耀斑觸發(fā)相關(guān)的臨界磁場特征.許多研究發(fā)現(xiàn),強(qiáng)磁場區(qū)域的極性反轉(zhuǎn)線在耀斑活動中起著重要作用.Sadykov等[1]分析了縱向磁場特征與太陽耀斑之間的關(guān)系,研究證實了極性反轉(zhuǎn)線活動區(qū)特征在耀斑預(yù)報過程中的獨(dú)特作用,并證明了僅使用縱向磁圖進(jìn)行耀斑預(yù)報的可能性.Cui等[2]從大量邁克爾遜多普勒成像儀(Michelson Doppler Imager,MDI)縱向磁圖中,計算了描述光球磁場特性的最大水平梯度、中性線長度和奇點(diǎn)數(shù)3個物理量,并擬合出了耀斑產(chǎn)生率與特征因子之間存在著S型函數(shù)的關(guān)系.他們的結(jié)果表明,太陽耀斑的產(chǎn)生率隨著磁場非勢性和復(fù)雜性的增加而增加.Yang等[3]的研究指出在耀斑之前沿著活動區(qū)的磁性中性線出現(xiàn)很強(qiáng)的剪切流動.Victor等[4]使用太陽磁場的極性反轉(zhuǎn)線(Polarity Inversion Line,PIL)長度、PIL附近的強(qiáng)磁場區(qū)域面積以及該區(qū)域的總通量作為太陽耀斑預(yù)測的特征,并獲得了良好的短期預(yù)測結(jié)果.Sharykin等[5]的觀測結(jié)果證明了位于色球?qū)拥入x子體中的主要能量釋放位置在具有強(qiáng)電流集中的PIL附近.開發(fā)用于檢測極性反轉(zhuǎn)線的方法為提高觀測數(shù)據(jù)使用效率、快速監(jiān)測太陽活動水平、改善空間天氣預(yù)報有著積極作用,同時也是一種提高我們對太陽上這些爆發(fā)現(xiàn)象認(rèn)識的方法.

極性反轉(zhuǎn)線位置通常由研究者根據(jù)研究經(jīng)驗人工繪制,或簡單使用等高線方法大致描繪.目前存在著幾類極性反轉(zhuǎn)線自動檢測方法,簡介如下:

基于像素空間卷積類方法.Volobuev等[6]使用一種元胞自動機(jī)的最近鄰檢測技術(shù),考慮磁圖的每個像素及其周圍的8個鄰域,若中心像素值與周圍鄰域像素值符號異號,則將其歸為極性反轉(zhuǎn)線上的像素.使用這個邏輯規(guī)則通過模板遍歷全圖找到所有極性反轉(zhuǎn)線集合.同樣,Steward等[7]將原始圖像逐個像素乘以從原始圖像向右移動一個像素獲得的新圖像,當(dāng)乘積結(jié)果發(fā)生異號并且大于一定的梯度閾值則歸為極性反轉(zhuǎn)線上的像素.對向上方向移動一個像素執(zhí)行相同的步驟,合并像素集獲得最終的極性反轉(zhuǎn)線.雖然此類方法較為簡單易于實現(xiàn),但從結(jié)果來看這些檢測的極性反轉(zhuǎn)線像素點(diǎn)較為零碎,受噪聲干擾較大.

基于梯度技術(shù)類檢測方法.Schrijver[8]為了識別強(qiáng)梯度的極性反轉(zhuǎn)線,首先對MDI磁圖構(gòu)造了二元正、負(fù)強(qiáng)磁場位圖,然后使用3×3的核進(jìn)行膨脹以創(chuàng)建膨脹的正負(fù)位圖,兩幅膨脹后的位圖的乘積不為零時被識別為強(qiáng)場區(qū)域的極性反轉(zhuǎn)線.Welsch等[9]使用Schrijver[8]的梯度檢測技術(shù)探索了關(guān)于極性反轉(zhuǎn)線附近強(qiáng)磁場梯度的起源,他們的結(jié)果表明新的磁場浮現(xiàn)確實可以導(dǎo)致強(qiáng)磁場梯度,但不是必要條件.

基于形態(tài)學(xué)擴(kuò)張技術(shù)類的檢測方法.Sadykov等[1]使用最早由Victor等[4]提出的磁圖分割的方法,將磁圖劃分為強(qiáng)正場、強(qiáng)負(fù)場和弱場(“中性”段)的區(qū)域.首先使用高斯濾波器對原始磁圖進(jìn)行平滑處理并應(yīng)用分割算法.然后,分別對正負(fù)區(qū)域應(yīng)用形態(tài)學(xué)擴(kuò)展程序,將正負(fù)片段擴(kuò)張后的交集作為極性反轉(zhuǎn)線集合.最后,我們過濾掉像素數(shù)少于一定數(shù)目的碎片.

這3類方法都存在著極性反轉(zhuǎn)線像素點(diǎn)較為細(xì)碎、不連續(xù)、受噪聲影響較大等問題.本文從支持向量機(jī)(Support Vector Machine,SVM)的模型出發(fā),將極性反轉(zhuǎn)線位置的探測問題轉(zhuǎn)化為模式識別中二分類問題.提出了一種基于支持向量機(jī)的極性反轉(zhuǎn)線檢測算法,自動探測與識別太陽動力學(xué)天文臺(Solar Dynamics Observatory,SDO)日震和磁成像儀(Helioseismic and Magnetic Imager,HMI)磁圖的極性反轉(zhuǎn)線位置.測試結(jié)果表明,此算法可以精確有效地檢測太陽活動區(qū)的極性反轉(zhuǎn)線.

2 算法原理

支持向量機(jī)是一種典型的二分類模型[10],主要用于模式識別中的分類問題.其基本思想是在特征空間中尋找一個最佳分割超平面,使得這個超平面能將樣本正確分類并且分類間距達(dá)到最大.支持向量機(jī)的一個關(guān)鍵概念就是引入了分類間隔,分類問題最終就可轉(zhuǎn)化成約束優(yōu)化問題求解.對一確定樣本,給定樣本點(diǎn)集合D={(x1,y1),(x2,y2),...,(xm,ym)},yi∈{?1,+1}.m為樣本點(diǎn)個數(shù).分類超平面可由下面表達(dá)式描述:

其中的w為超平面的法向量,描述了分割超平面的方向信息,b為位移,x為樣本點(diǎn).超平面由(w,b)唯一決定.樣本空間中的任一點(diǎn)x到超平面的距離可表示為:

有了超平面和分類間隔的描述就可以確定一個分類器了.假設(shè)分割超平面能將樣本正確劃分,對于(xi,yi)∈D,二分類器描述為:

即:

距離超平面最近的點(diǎn)使得等號成立,叫做“支持向量”.兩個異側(cè)支持向量到超平面距離之和稱為間隔γ,用下式表示:

要找一個劃分超平面能將樣本正確分類并且分類間隔達(dá)到最大,也就是找到(w,b)使其(5)式分類器成立,即:

支持向量機(jī)基本型便可寫為:

在實際應(yīng)用中還存在錯誤分類的問題,僅通過最大間距找到?jīng)Q策邊界是不夠的,還需要考慮錯誤分類的情況.因此需要一個懲罰因子C,增強(qiáng)模型的容錯能力.它在優(yōu)化函數(shù)里主要是平衡支持向量的復(fù)雜度和誤分類率這兩者之間的關(guān)系.

實際問題中,我們的數(shù)據(jù)在原始樣本空間中并不一定線性可分,無法找到這樣一個線性超平面將樣本正確分類.對于線性不可分問題,SVM使用一種叫做“核函數(shù)”的技術(shù),將原始線性不可分低維空間映射到一個更高維度特征空間使其線性可分.

設(shè)φ(x i)為映射后的特征向量,則非線性優(yōu)化問題轉(zhuǎn)化為:

其中φ(x i)為低維到高維的映射核函數(shù).常用核函數(shù)見表1.

表1 常用核函數(shù)及其表達(dá)式Table 1 Com m on ker nel functions and their ex p ressions

表中,x j為D中樣本點(diǎn),d為多項式次數(shù),σ是高斯核函數(shù)帶寬,tanh為雙曲正切函數(shù).β、θ為根據(jù)實驗人工設(shè)置的參數(shù),β是分類類別的倒數(shù),θ是偏移量.

本文提出的基于支持向量機(jī)的極性反轉(zhuǎn)線檢測算法,關(guān)鍵在于如何轉(zhuǎn)換成二分類問題,下面將詳細(xì)描述具體算法流程:

步驟1:待檢測區(qū)域截取,活動區(qū)的截取直接使用鼠標(biāo)點(diǎn)選的方式;

步驟2:對截取后的圖像進(jìn)行預(yù)處理,包括平滑去噪等;

步驟4:SVM輸出,將每個像素點(diǎn)值的正負(fù)代數(shù)符號轉(zhuǎn)換成輸出標(biāo)簽,用+1或者?1表示;

步驟5:核函數(shù)映射,選取徑向基函數(shù)(Radial Basis Function,RBF)作為核函數(shù),映射到特征空間;

步驟6:超平面計算,通過選取的特征向量與標(biāo)簽利用SVM計算(1)式中的w與b;

步驟7:極性反轉(zhuǎn)線位置確定,有了分割超平面參數(shù)w與b即確定了極性反轉(zhuǎn)線位置.

本文算法的具體流程如圖1所示.

3 實驗過程

本文的實驗數(shù)據(jù)來自于SDO/HMI,數(shù)據(jù)形式為FITS格式.圖2所示為2017年9月5日NOAA12673活動區(qū)縱向磁圖.我們以圖像左上點(diǎn)為原點(diǎn),從原點(diǎn)自左向右為X軸,自上向下為Y軸建立坐標(biāo)系,圖像分辨率大小為160像素×160像素.

圖1 基于支持向量機(jī)的極性反轉(zhuǎn)線檢測算法流程Fig.1 Flow chart of p olarity inversion line detection algorithm based on supp ort vector machine

基于支持向量機(jī)的極性反轉(zhuǎn)線檢測算法具體實驗流程如下:

(1)磁場活動區(qū)區(qū)域截取、圖像去噪、平滑等預(yù)處理.

會上，中國農(nóng)業(yè)大學(xué)教授、資源與環(huán)境學(xué)院副院長江榮風(fēng)，安徽農(nóng)業(yè)大學(xué)黨委常委、副校長姚佐文，安徽六國化工副總經(jīng)理、總工程師沈浩，代表三方簽署戰(zhàn)略合作協(xié)議。

從全日面的太陽磁場圖像中截取所需要檢測的局部活動區(qū),然后,使用高斯濾波器對磁場圖像進(jìn)行平滑去噪處理,該操作能去除很多無關(guān)細(xì)碎的小塊,保持極性反轉(zhuǎn)線的連續(xù)性,增強(qiáng)檢測結(jié)果的準(zhǔn)確性與穩(wěn)定性.

(2)二分類模型的轉(zhuǎn)換.

將圖像數(shù)據(jù)轉(zhuǎn)換成可以輸入到支持向量機(jī)模型中的向量.基于極性反轉(zhuǎn)線的概念,極性反轉(zhuǎn)線的位置只與正負(fù)磁極位置相關(guān).具體做法是,遍歷每個像素點(diǎn),記錄其圖像坐標(biāo)p(u,v)作為SVM特征向量輸入.將每個像素點(diǎn)灰度值的正負(fù)代數(shù)符號轉(zhuǎn)換成輸出類別,用+1或者?1表示.因此輸入向量為(u,v,±1).我們將大于一定閾值?的像素歸為正極性類,用+1表示,而將小于閾值?的像素歸為負(fù)極性類,用?1表示,這樣就將轉(zhuǎn)換成一個二分類問題.這里閾值?的選擇很重要,因為小的閾值將導(dǎo)致PIL的細(xì)碎離散片段增加,這些小的片段,對研究者并無實際用途.較大的閾值將不能正確完全檢測PIL.在本實驗中圖像分辨率大小為160像素×160像素,通過大量實驗,結(jié)果表明30–50 Gs能夠取得較好效果.

(3)構(gòu)建支持向量機(jī)模型.

模型構(gòu)建根據(jù)支持向量機(jī)模型,定義規(guī)則,正極性類yi=+1,對所有正極性樣本點(diǎn)以(u,v)表示,使得:

負(fù)極性類yi=?1,對所有負(fù)極性樣本點(diǎn)以(u,v)表示,使得:

由于我們的數(shù)據(jù)線性不可分,需要選擇一個合適的核函數(shù).對于非線性可分的問題,高斯核函數(shù)也稱為RBF核函數(shù)能夠取得較好結(jié)果.高斯核如下式:

其中,高斯核σ參數(shù)取值與樣本的劃分精細(xì)程度有關(guān):σ越大,低維空間中選擇的曲線越復(fù)雜,分的類別越細(xì),容易出現(xiàn)過擬合.σ越小,分的類別越粗,可能導(dǎo)致無法將數(shù)據(jù)區(qū)分開來,容易出現(xiàn)欠擬合.

給定C,增強(qiáng)模型的容錯能力.它在優(yōu)化函數(shù)里主要是平衡支持向量的復(fù)雜度和誤分類率這兩者之間的關(guān)系.當(dāng)C比較大時,我們的損失函數(shù)也會越大,這樣我們會有更多的支持向量,也就是說支持向量和超平面的模型也會變得越復(fù)雜,也容易過擬合.反之,當(dāng)C比較小時,會選擇較少的樣本來做支持向量,最終的支持向量和超平面的模型也會簡單.

(4)計算分割超平面,將決策邊界在磁圖中畫出.

將第3步轉(zhuǎn)換后的數(shù)據(jù)輸入到支持向量機(jī)模型中,計算分割超平面求解(w,b)使(wTx+b)=1,將計算出的決策邊界在磁圖上畫出.如圖3所示是使用高斯核函數(shù)參數(shù)C=0.1,σ=0.05最終檢測結(jié)果.

圖2 2017年9月5日NOAA AR(Active Region)12673磁圖Fig.2 Magnetograms of NOAA AR12673 on 2017 September 5

圖3 NOAA AR12673磁圖,綠色的線為閾值為30 Gs,使用高斯核函數(shù)參數(shù)C=0.1,σ=0.05極性反轉(zhuǎn)線檢測結(jié)果.Fig.3 Magnetograms of NOAA AR12673,the green line is the p olarity inversion lines detection result using a threshold of 30 Gs and a Gaussian kernel function parameter of C=0.1,σ=0.05.

根據(jù)研究者的目的不同,可以使用不同的參數(shù)控制檢測的精細(xì)程度,如圖4所示.

為了說明算法的有效性,我們對比基于像素空間卷積類方法,此方法也是較為常見的方法.對同樣的數(shù)據(jù),使用模板計算每個像素與其周圍的4個鄰域像素的卷積,若中心像素值與周圍領(lǐng)域像素值符號異號,則將其歸為極性反轉(zhuǎn)線上的像素.使用這個邏輯規(guī)則通過模板遍歷全圖找到所有極性反轉(zhuǎn)線集合,如圖5所示為使用閾值為50 Gs的實驗結(jié)果.

從檢測效果來看,此方法雖然簡單易于實現(xiàn),但極性反轉(zhuǎn)線像素點(diǎn)較為細(xì)碎,不連續(xù),受噪聲影響較大.使用不同的閾值時,檢測結(jié)果會存在很大差異,例如使用稍小的閾值就會引入大量寧靜區(qū)的污染像素,造成極性反轉(zhuǎn)線統(tǒng)計結(jié)果有較大差異.同時對于不同的磁圖,閾值選擇無法統(tǒng)一標(biāo)準(zhǔn),人為干預(yù)性較大.

極性反轉(zhuǎn)線的長度通常是研究者所關(guān)心的問題,本文算法在檢測結(jié)果上能較為直觀地反映極性反轉(zhuǎn)線的真實位置,同時極性反轉(zhuǎn)線的連續(xù)性好.極性反轉(zhuǎn)線長度可定義為極性反轉(zhuǎn)線所占像素數(shù)之和.因此在后續(xù)計算極性反轉(zhuǎn)線長度時可以在檢測結(jié)果的基礎(chǔ)上采用如上定義統(tǒng)計最長的極性反轉(zhuǎn)線像素之和.

圖4 NOAA AR12673磁圖,綠色的線為使用閾值為50 Gs,高斯核函數(shù)參數(shù)為C=0.1,σ=0.1的極性反轉(zhuǎn)線檢測結(jié)果.Fig.4 Magnetograms of NOAA AR12673,the green line is the polarity inversion lines detection result using a threshold of 50 Gs and a Gaussian kernel function parameter of C=0.1,σ=0.1.

圖5 NOAA AR12673磁圖,綠色點(diǎn)為極性反轉(zhuǎn)線上的像素.Fig.5 Magnetograms of NOAA AR12673,green colors show the pixels of p olarity inversion lines.

4 總結(jié)與討論

本文從支持向量機(jī)的模型出發(fā),結(jié)合極性反轉(zhuǎn)線的特性,將極性反轉(zhuǎn)線位置的探測問題轉(zhuǎn)化為一個模式識別中的二分類問題.提出了一種基于支持向量機(jī)的有效算法,以太陽動力學(xué)觀測臺日震及磁成像儀磁圖為實驗數(shù)據(jù)檢測極性反轉(zhuǎn)線.實驗結(jié)果表明,此算法可以有效檢測太陽活動區(qū)的極性反轉(zhuǎn)線,檢測精度較高,檢測結(jié)果也可疊加在同一區(qū)域不同波段的太陽圖像上,便于研究者開展后續(xù)工作.另外,基于支持向量機(jī)原理,此方法也可推廣至寧靜區(qū)大尺度暗條極性反轉(zhuǎn)線的檢測.今后將會進(jìn)一步優(yōu)化算法應(yīng)對全局區(qū)域的檢測.