蔣 博 劉 磊 鄭 勝 楊珊珊 曾曙光 黃 瑤 羅驍域

(三峽大學(xué)天文與空間科學(xué)研究中心 宜昌 443002)

1 引言

太陽表面存在各種不同的活動(dòng)現(xiàn)象, 它們與太陽大氣中磁場的分布和位形密切相關(guān), 涉及等離子體的加熱和運(yùn)動(dòng)、電磁輻射的增強(qiáng)、高能粒子的加速和波動(dòng)等各種物理過程[1]. 太陽黑子是太陽表面最基本的活動(dòng)現(xiàn)象, 黑子的多少表征著太陽活動(dòng)程度的高低, 由黑子組成的活動(dòng)區(qū)是太陽爆發(fā)(太陽耀斑和日冕物質(zhì)拋射)的源區(qū). 它驅(qū)動(dòng)著太陽活動(dòng), 對近地空間、行星際和日球環(huán)境有著極為重要的影響[2]. 因此, 對于太陽活動(dòng)區(qū)域進(jìn)行準(zhǔn)確的識別是一項(xiàng)重要的工作[3].

太陽活動(dòng)區(qū)的檢測是一種典型的多目標(biāo)檢測問題,已有很多學(xué)者對太陽活動(dòng)區(qū)的檢測進(jìn)行研究.早期的檢測通常采用經(jīng)典的圖像處理技術(shù), 例如閾值分割法[4]、區(qū)域生長法[5–6]等, 但是這些方法的檢測效果強(qiáng)烈依賴于強(qiáng)度閾值和區(qū)域生長的邊界閾值等參數(shù)的設(shè)置[7]. Zhang等[8]提出的自動(dòng)圖像處理與機(jī)器學(xué)習(xí)相結(jié)合的自動(dòng)檢測算法會因?yàn)榛顒?dòng)區(qū)的重疊導(dǎo)致無法正確聚類. Caballero等[9]提出的與極紫外成像望遠(yuǎn)鏡和日光層天文臺航天器相結(jié)合的自動(dòng)檢測系統(tǒng)以及Higgins等[5]提出的使用SOHO (Solar and Heliospheric Observatory)邁克爾遜多普勒成像儀(Michelson Doppler Imager,MDI)與太陽監(jiān)測器活動(dòng)區(qū)域跟蹤(The Solar Monitor Active Region Tracking, SMART)算法相結(jié)合的自動(dòng)檢測系統(tǒng)均容易導(dǎo)致多個(gè)相鄰的太陽活動(dòng)區(qū)被誤檢測為一個(gè)或一個(gè)太陽活動(dòng)區(qū)被誤檢測為多個(gè). 近年來, 人工智能在天文學(xué)上的應(yīng)用與日俱增. 例如, 朱健等[10]基于YOLOv3 (You Only Look Once)和DeepSort算法對太陽活動(dòng)區(qū)進(jìn)行的檢測和跟蹤, 采用16361個(gè)太陽活動(dòng)區(qū)樣本進(jìn)行訓(xùn)練, 達(dá)到了92%的檢測準(zhǔn)確率. LeNet-5以及AlexNet等卷積神經(jīng)網(wǎng)絡(luò)也被廣泛地用于天文學(xué)研究中[11–12]. 這些人工智能方法在天文學(xué)研究中雖然有著較高的準(zhǔn)確率, 但是需要提前對大量樣本數(shù)據(jù)進(jìn)行訓(xùn)練與測試后才能使用.

尺度不變特征變換(Scale-Invariant Feature Transform, SIFT)是一種基于局部描述子的特征點(diǎn)檢測方法, 對于光線、噪聲、視角改變的容忍度相當(dāng)高[13–14], 對于含有許多內(nèi)部結(jié)構(gòu)、特征點(diǎn)較明顯的太陽磁場圖像具有較好的檢測效果.2018年,楊盼等[15]采用SIFT算法將懷柔太陽觀測基地的局部太陽磁場圖像與日震和磁像儀(Helioseismic and Magnetic Imager, HMI)的全日面太陽磁場圖像進(jìn)行配準(zhǔn)與定位. 近年來, 聚類算法被廣泛地應(yīng)用在目標(biāo)識別中. 其中, 密度峰值聚類(Clustering by Fast Search and Find of Density Peaks, DPC)算法是一種基于數(shù)據(jù)點(diǎn)密度進(jìn)行自動(dòng)分類的聚類算法[16]. 該算法基于聚類中心擁有較高的密度及聚類中心點(diǎn)之間距離較大兩個(gè)假設(shè), 能夠直觀地得到聚類的個(gè)數(shù)、自動(dòng)地發(fā)現(xiàn)聚類中心并且能排除異常點(diǎn)同時(shí)不受聚類對象嵌入空間維數(shù)的限制, 實(shí)現(xiàn)任意形狀數(shù)據(jù)的高效聚類.

本文針對太陽活動(dòng)區(qū)結(jié)構(gòu)較為復(fù)雜的特點(diǎn), 采用SIFT特征提取算法提取太陽活動(dòng)區(qū)的特征點(diǎn),利用密度峰值聚類算法對提取的特征點(diǎn)進(jìn)行聚類,實(shí)現(xiàn)了太陽活動(dòng)區(qū)的自動(dòng)檢測. 該算法可以在不需要人工交互的情況下有效檢測并識別太陽活動(dòng)區(qū).

2 檢測算法原理

SIFT和DPC算法主要原理如下.

2.1 SIFT算法

SIFT是一種計(jì)算機(jī)視覺的算法. 它用來偵測與描述影像中的局部性特征, 它在空間尺度中尋找極值點(diǎn), 并提取出其位置、尺度、旋轉(zhuǎn)不變量, 此算法由Lowe在1999 年發(fā)表[14],2004年完善總結(jié)[13].該方法步驟如下:

(1)尺度空間的生成

尺度空間的理論方法是檢測特征不變性的主要基礎(chǔ). 尺度空間理論的目的是模擬圖像的多尺度特征, Koenderink[17]證明了高斯核函數(shù)是實(shí)現(xiàn)尺度變換的唯一變換核. 二維圖像與高斯核函數(shù)的卷積可以得到該圖像在不同尺度下的尺度空間L(x,y,σ).

式中,x、y為像素坐標(biāo),I(x,y)為灰度值,G(x,y,σ)為高斯核函數(shù),σ是高斯正態(tài)分布的方差. 高斯核函數(shù)的定義如下:

其中, 正態(tài)分布方差σ在圖像中叫做尺度空間因子,它反映了圖像被平滑的程度, 尺度越小表示圖像被平滑得越小, 對應(yīng)于圖像的細(xì)節(jié)特征, 大尺度對應(yīng)于圖像的概貌特征.

(2)檢測尺度空間的極值點(diǎn)

為高效地在尺度空間檢測出穩(wěn)定的極值,Lowe[13]使用尺度空間中差分高斯DOG的極值作為判斷依據(jù). DOG算子定義如下:

其中,k是相鄰兩個(gè)尺度空間的比例因子.

為了確保特征點(diǎn)的穩(wěn)定性和匹配準(zhǔn)確性, 對SIFT候選點(diǎn)集進(jìn)行篩選, 去除其中對比度低的特征點(diǎn)與圖像邊緣的特征點(diǎn), 增強(qiáng)匹配的穩(wěn)定性. 以上兩步驟的示意圖如圖1所示.

圖1 尺度變換與極值點(diǎn)檢測[14]Fig.1 Scale transformation and extreme detection[14]

(3)計(jì)算方向幅值

為了區(qū)分不同關(guān)鍵點(diǎn)的屬性, 可以利用其鄰域像素的梯度變化特性來計(jì)算關(guān)鍵點(diǎn)附近各個(gè)方向的方向幅值, 并繪制該關(guān)鍵點(diǎn)附近的方向累積直方圖, 其中幅值累積最大的方向?yàn)樘卣鼽c(diǎn)的主梯度方向.

(4)關(guān)鍵點(diǎn)描述

以關(guān)鍵點(diǎn)周圍8×8的窗口為例, 將坐標(biāo)軸旋轉(zhuǎn)為主梯度方向, 計(jì)算關(guān)鍵點(diǎn)周圍8×8的窗口中每一個(gè)像素的梯度, 并且使用高斯下降函數(shù)降低遠(yuǎn)離中心的關(guān)鍵點(diǎn)的權(quán)重. 這樣就能生成一個(gè)2×2×8= 32維的描述子, 每一維都可以表示2×2格子中的一個(gè)梯度的方向與尺度大小特征, 如圖2所示.

圖2 生成關(guān)鍵點(diǎn)描述子示意圖[14]Fig.2 Generate the keypoint descriptor[14]

2.2 密度峰值聚類

基于DPC算法可以識別任意大小、形狀的聚類對象, 同時(shí)可以有效地濾除噪聲[16]. DPC算法2014年發(fā)表在Science雜志上, 使得基于密度的聚類方法取得質(zhì)的飛躍. DPC的核心思想是認(rèn)為聚類中心的密度顯著大于其鄰域內(nèi)其他點(diǎn)的密度, 且聚類中心之間的距離相對較遠(yuǎn). DPC的有效檢測建立在兩個(gè)假設(shè)上: (1)聚類中心被局部密度較低的點(diǎn)包圍; (2)聚類中心點(diǎn)與其他聚類中心點(diǎn)的距離相對較大. DPC算法的步驟如下.

(1)使用截?cái)嗪?Cut offkernel)計(jì)算每個(gè)點(diǎn)的密度ρi:

(2)根據(jù)每個(gè)數(shù)據(jù)點(diǎn)周圍一定范圍內(nèi)數(shù)據(jù)點(diǎn)的數(shù)量, 計(jì)算得到數(shù)據(jù)點(diǎn)的密度{ρ1,ρ2,··· ,ρn},n為點(diǎn)數(shù), 并對其進(jìn)行排序, 計(jì)算每個(gè)點(diǎn)的距離δi:

當(dāng)存在密度高于當(dāng)前數(shù)據(jù)點(diǎn)的數(shù)據(jù)點(diǎn)時(shí), 當(dāng)前數(shù)據(jù)點(diǎn)的距離為比該點(diǎn)密度大的所有數(shù)據(jù)點(diǎn)距離的最小值. 其中, 密度最大的數(shù)據(jù)點(diǎn)的距離指定為所有數(shù)據(jù)點(diǎn)中距離最大值.

(3)通過設(shè)定密度和距離的閾值, 選取聚類中心點(diǎn). 根據(jù)非類中心和類中心之間的距離, 將其他非類中心點(diǎn)分配到距離最近的聚類中心點(diǎn), 從而完成聚類.

3 算法設(shè)計(jì)及實(shí)驗(yàn)分析

本文采用的觀測數(shù)據(jù)來自于太陽動(dòng)力學(xué)天文臺(Solar Dynamics Observatory,SDO)上HMI的全日面縱向磁圖[17–18]. 實(shí)驗(yàn)使用的所有全日面圖像數(shù)據(jù)從http://jsoc.stanford.edu/ajax/lookdata.html網(wǎng)站上以fits格式下載, 將每張樣本圖像轉(zhuǎn)換成2048×2048像素的JPG圖像, 使用SIFT和DPC算法相結(jié)合的方法來實(shí)現(xiàn)太陽活動(dòng)區(qū)的自動(dòng)檢測, 實(shí)驗(yàn)流程及結(jié)果如下.

3.1 算法設(shè)計(jì)

基于SIFT與DPC的太陽活動(dòng)區(qū)自動(dòng)檢測算法流程圖如圖3所示, 具體步驟如下:

圖3 太陽活動(dòng)區(qū)自動(dòng)檢測算法流程Fig.3 Flow chart of solar activity regions automatic detection

(1)圖像預(yù)處理, 對原始太陽磁場圖像進(jìn)行對比度增強(qiáng), 增強(qiáng)活動(dòng)區(qū)的特征信息;

(2)基于SIFT算法, 對增強(qiáng)后的圖像進(jìn)行特征點(diǎn)提取, 獲得該全日面特征點(diǎn)集, 根據(jù)拍攝時(shí)間, 將特征點(diǎn)的像素坐標(biāo)轉(zhuǎn)化為經(jīng)緯度坐標(biāo);

(3)對提取出來的特征點(diǎn)進(jìn)行密度峰值聚類,根據(jù)特征點(diǎn)的密度分布與特征點(diǎn)相互間的距離特征, 可自動(dòng)將特征點(diǎn)聚類;

(4)根據(jù)聚類決策圖, 采用一定的判據(jù)(經(jīng)統(tǒng)計(jì),聚類中心點(diǎn)的密度需大于所有特征點(diǎn)密度的均值減0.5倍方差,距離需大于所有特征點(diǎn)距離的均值加3倍方差, 且該特征點(diǎn)在決策圖上到坐標(biāo)原點(diǎn)的距離需大于決策圖對角線長度的一半)選出候選區(qū)域,完成太陽活動(dòng)區(qū)的初步定位;

(5)基于最大類間方差法[19]對初步檢測結(jié)果進(jìn)行自動(dòng)圖像分割, 再用面積濾波去除噪聲之后, 獲取所得活動(dòng)區(qū)的最小外接矩形, 即可完成對太陽活動(dòng)區(qū)的精確定位;

(6)對日面經(jīng)度延伸較大(例如大于20°)的活動(dòng)區(qū)進(jìn)行二次檢測, 根據(jù)二次檢測決策圖, 采用二次聚類的判據(jù)(經(jīng)統(tǒng)計(jì), 聚類中心點(diǎn)的密度需大于所有特征點(diǎn)密度的均值加方差, 距離需大于所有特征點(diǎn)距離的均值加6倍方差)對其進(jìn)行判別. 此時(shí), 如果滿足上述判據(jù)的特征點(diǎn)數(shù)量為2, 則需根據(jù)磁通量進(jìn)一步判斷, 若兩特征點(diǎn)對應(yīng)區(qū)域中均只含有一種極性, 且這兩區(qū)域的極性相反, 則它們?yōu)橥换顒?dòng)區(qū); 反之, 則應(yīng)劃為兩個(gè)活動(dòng)區(qū).

3.2 實(shí)驗(yàn)結(jié)果與分析

根據(jù)上述算法思想, 我們編寫了太陽活動(dòng)區(qū)檢測程序. 實(shí)驗(yàn)結(jié)果如下, 圖4 (a)是2017年1月24日的全日面磁圖, 圖4 (b)顯示的是對其進(jìn)行對比度增強(qiáng)后的結(jié)果.

圖4 原始太陽圖像與預(yù)處理結(jié)果. (a)原始圖像; (b)預(yù)處理結(jié)果圖像.Fig.4 The original image and image after preprocessing. (a)The original image; (b) image after preprocessing.

然后, 使用SIFT算法提取增強(qiáng)后磁場圖像的特征點(diǎn), 結(jié)果如圖5所示. 在這張圖中我們總共檢測到了1427個(gè)特征點(diǎn), 可以看出在太陽活動(dòng)區(qū)周圍的特征點(diǎn)分布較為集中, 顯著多于寧靜區(qū)的特征點(diǎn),檢測結(jié)果很好地代表了太陽活動(dòng)區(qū)的特征, 為下一步利用聚類算法進(jìn)行活動(dòng)區(qū)檢測提供了基礎(chǔ).

圖5 SIFT算法檢測到的特征點(diǎn)分布Fig.5 The distribution of feature points detected by the SIFT algorithm

圖6 展示了上述1427個(gè)特征點(diǎn)的密度峰值聚類結(jié)果圖. 圖6 (a)為密度峰值聚類決策圖, 其橫坐標(biāo)ρ為特征點(diǎn)歸一化后的密度, 縱坐標(biāo)δ表示特征點(diǎn)和離它最近且密度大于該點(diǎn)的特征點(diǎn)之間歸一化后的距離. 因此特征點(diǎn)在決策圖中的位置越靠近右上, 表示該點(diǎn)越靠近聚類中心, 反之亦然. 通過統(tǒng)計(jì)分析, 我們發(fā)現(xiàn)太陽活動(dòng)區(qū)的聚類中心應(yīng)滿足條件: (1)密度大于所有特征點(diǎn)的密度均值減0.5倍方差; (2)距離大于所有特征點(diǎn)的距離均值加3倍方差;(3)距離坐標(biāo)原點(diǎn)的距離大于對角線長度的一半,即. 在圖6 (a)的DPC決策圖中, 所有特征點(diǎn)的密度均值為0.26, 密度方差為0.23, 距離均值為0.03,距離方差為0.07; 相應(yīng)的密度判據(jù)為0.145, 距離判據(jù)為0.24, 距離原點(diǎn)的距離判據(jù)為滿足上述判據(jù)的特征點(diǎn)一共有兩個(gè), 在圖6 (a)中分別用藍(lán)色和青色方框加以了標(biāo)注. 由此可以判斷出該磁場圖像中的特征點(diǎn)有兩個(gè)聚類中心. 圖6 (b)為特征點(diǎn)的聚類結(jié)果圖, 其中Dim1與Dim2分別表示圖像的橫縱坐標(biāo)(單位為像素). 在該圖中, 聚類算法將該圖像的特征點(diǎn)自動(dòng)分為了兩類, 聚類中心分別用深藍(lán)色和淺藍(lán)色菱形加以標(biāo)注; 右邊紅色點(diǎn)區(qū)域?yàn)榈?類, 左邊黑色區(qū)域?yàn)榈?類, 綠色點(diǎn)區(qū)域?yàn)樵肼?

根據(jù)該聚類結(jié)果在全日面磁圖中的像素位置,可實(shí)現(xiàn)對太陽活動(dòng)區(qū)在全日面磁圖中的初步定位.對初步定位的結(jié)果, 使用最大類間方差[19]進(jìn)行圖像自動(dòng)分割, 并使用面積濾波去除噪聲, 即可獲得太陽活動(dòng)區(qū)的精確定位. 對全日面磁圖的太陽活動(dòng)區(qū)檢測結(jié)果如圖7所示, 其中藍(lán)色虛線框是對太陽活動(dòng)區(qū)初步定位的檢測結(jié)果, 紅色實(shí)線框是對太陽活動(dòng)區(qū)精確檢測的結(jié)果.

圖7 密度峰值聚類算法檢測結(jié)果Fig.7 The results of the proposed algorithm

為驗(yàn)證算法的有效性, 我們對不同時(shí)間拍攝的全日面太陽磁圖進(jìn)行了檢測, 圖8是部分太陽磁圖的檢測結(jié)果, 其中, 圖8 (a)日面上僅存在一個(gè)太陽活動(dòng)區(qū), 圖8 (b)日面上有3個(gè)太陽活動(dòng)區(qū). 結(jié)果表明, 無論全日面太陽磁圖中存在單個(gè)或是多個(gè)太陽活動(dòng)區(qū), 該算法均能夠有效地檢測出太陽活動(dòng)區(qū).

圖8 不同時(shí)間太陽磁場圖像的自動(dòng)檢測結(jié)果Fig.8 The results of automatic detection for different solar magnetic field images

當(dāng)兩個(gè)活動(dòng)區(qū)位置相近時(shí), 有可能會產(chǎn)生誤檢, 即相近區(qū)域的兩個(gè)活動(dòng)區(qū)被誤檢為一個(gè)活動(dòng)區(qū). 為進(jìn)一步提高檢測準(zhǔn)確性, 我們提出了對誤檢區(qū)域進(jìn)行二次檢測的方法. 對精確定位的結(jié)果進(jìn)行檢測, 若日面經(jīng)度延伸大于20°, 則判定該區(qū)域可能存在兩個(gè)活動(dòng)區(qū). 通過對提取出的誤檢區(qū)域使用SIFT算法與密度峰值聚類的二次檢測, 可實(shí)現(xiàn)對誤檢區(qū)域的修正, 檢測出正確的太陽活動(dòng)區(qū). 從而有效解決將多個(gè)活動(dòng)區(qū)誤檢為一個(gè)活動(dòng)區(qū)的問題,提高檢測的準(zhǔn)確性.

圖9展示了本文方法與區(qū)域生長法的對比結(jié)果圖. 圖9 (a)是本文方法的檢測結(jié)果, 其中藍(lán)色虛線框是第1次檢測結(jié)果, 紅色實(shí)線框?yàn)槎螜z測的結(jié)果, 圖9 (b)顯示了對于該圖像基于區(qū)域生長法的太陽活動(dòng)區(qū)檢測結(jié)果圖. 由于圖9 (a)中左上角區(qū)域的精確檢測定位結(jié)果日面經(jīng)度延伸大于20°, 因此對該區(qū)域進(jìn)行二次檢測. 圖10 (a)是該區(qū)域的二次檢測聚類決策圖, 在該圖中, 滿足二次聚類判據(jù)的特征點(diǎn)共有兩個(gè), 分別用編號1、2加以標(biāo)注; 圖10(b)是根據(jù)決策圖的決策結(jié)果對上述區(qū)域進(jìn)行二次檢測,得到的兩個(gè)活動(dòng)區(qū)分別用編號1、2加以標(biāo)注.這樣一來, 一次檢測中多個(gè)鄰近活動(dòng)區(qū)被誤檢為一個(gè)活動(dòng)區(qū)的問題得到了解決.

圖9 (a)本文算法檢測結(jié)果; (b)區(qū)域生長法檢測結(jié)果.Fig.9 (a) The results of the proposed algorithm; (b) the results of the region growing algorithm.

圖10 針對誤檢的太陽活動(dòng)區(qū)二次檢測結(jié)果. (a)密度峰值聚類決策圖; (b)精確定位區(qū)域的二次檢測結(jié)果.Fig.10 The second detection results of erroneous solar activity regions. (a) Decision graph of second detection; (b) the second detection results of accurate positioning.

通過兩種方法的對比, 本文的檢測算法具有以下4個(gè)優(yōu)勢: (1)自動(dòng)化程度高, 相比區(qū)域生長法,本文算法不需要提前選取種子點(diǎn), 所有特征點(diǎn)的提取全部是自動(dòng)進(jìn)行的; (2)區(qū)域生長法可能將相鄰的活動(dòng)區(qū)檢測為一個(gè)活動(dòng)區(qū), 比如圖9 (a)中編號為1和2的兩塊活動(dòng)區(qū), 在圖9 (b)中區(qū)域生長法誤檢測成了一個(gè)活動(dòng)區(qū)(編號1); (3)對于較大型活動(dòng)區(qū)的延伸區(qū)域以及雙極相隔較遠(yuǎn)的太陽活動(dòng)區(qū), 區(qū)域生長法可能會導(dǎo)致一個(gè)活動(dòng)區(qū)被誤檢為兩個(gè), 比如圖9 (a)中編號3和4的兩塊活動(dòng)區(qū), 在圖9(b)中區(qū)域生長法分別將其誤檢成了編號2和3以及編號4和5兩個(gè)活動(dòng)區(qū); (4)區(qū)域生長法的活動(dòng)區(qū)候選區(qū)域通過強(qiáng)度閾值確定, 因此會對強(qiáng)磁性但非活動(dòng)區(qū)的區(qū)域產(chǎn)生誤檢[10].

經(jīng)過上述結(jié)果對比分析, 本文算法可以有效地自動(dòng)檢測全日面上的太陽活動(dòng)區(qū). 而傳統(tǒng)方法中存在的由于活動(dòng)區(qū)雙極相距較遠(yuǎn)而導(dǎo)致被誤檢為多個(gè)活動(dòng)區(qū), 或者多個(gè)相距較近的活動(dòng)區(qū)被誤檢測為一個(gè)活動(dòng)區(qū)等問題[10], 本文算法也能較好地解決.

4 總結(jié)與展望

本文提出了一種基于SIFT特征檢測與密度峰值聚類相結(jié)合的太陽活動(dòng)區(qū)自動(dòng)檢測算法. 該算法分為兩步: 首先, 采用SIFT算法提取全日面圖像的特征點(diǎn); 其次, 基于密度峰值聚類算法對特征點(diǎn)進(jìn)行聚類, 初步確定活動(dòng)區(qū), 然后采用自動(dòng)閾值分割對活動(dòng)區(qū)精確定位. 對不同復(fù)雜程度的全日面磁場圖像進(jìn)行了活動(dòng)區(qū)檢測, 結(jié)果表明該算法可以較為準(zhǔn)確地提取出太陽活動(dòng)區(qū), 且魯棒性較強(qiáng), 自動(dòng)化程度較高. 需要指出的是, 該算法也有不足之處, 對太陽邊緣的活動(dòng)區(qū)進(jìn)行檢測時(shí)的準(zhǔn)確率會下降, 因此如何提高邊緣區(qū)域的活動(dòng)區(qū)檢測準(zhǔn)確度將是我們下一步研究的方向.