摘要：日冕物質拋射(CoronalMassEjection,CME)的檢測是建立CME事件庫和實現對CME在行星際傳播的預報的重要前提.透過VisualGeometryGroup(VGG)16卷積神經網路方法對日冕儀影象進行自動分類.基於大角度光譜日冕儀C2的白光日冕儀影象,根據是否觀測到CME對影象進行標記.將標記分類的資料集用於VGG模型的訓練,該模型在測試集分類的準確率達到92.5%.根據檢測得到的標籤結果,結合時空連續性規則,消除了誤判區域,有效分類出CME影象序列.與CoordinatedDataAnalysisWorkshops(CDAW)人工事件庫比較,分類出的CME影象序列能夠較完整地包含CME事件,且對弱CME結構有較高的檢測靈敏度.未來先進天基太陽天文臺(AdvancedSpace-basedSolarObservatory,ASO-S)衛星的萊曼阿爾法太陽望遠鏡將搭載有白光日冕儀(SolarCoronaImager,SCI),使用此分類方法將該儀器產生的日冕影象按有無CME分類.含CME標籤的影象將推送給中國的各空間天氣預報中心,對CME進行預警.
關鍵詞： 影象處理 太陽 技術 資料分析 日冕物質拋射
1、引言
日冕物質拋射(CoronalMassEjection,CME)是太陽大氣中最劇烈、尺度最大的活動現象,表現為在短時間內日冕結構發生明顯的變化,並伴有1011–1013kg攜帶磁場的等離子體拋射.當日冕物質拋射的方向朝著地球時,可能會與地球磁層發生相互作用,引起近地空間的地磁暴、極光等現象,會對通訊系統和電力系統等產生干擾,嚴重時會造成巨大的經濟損失.因此,CME到達地球的實時預報對空間天氣環境的監測十分重要.
CME的自動標註和檢測是實現CME預報的重要前提.太陽和日球層天文臺(SolarandHeliosphericObservatory,SOHO)搭載的大角度光譜日冕儀(LargeAngleandSpectrometricCoronagraphExperiment,LASCO)能夠觀測太陽日冕活動.LASCO由3臺視場不同的日冕儀構成,其中LASCOC2視場的範圍大約是太陽直徑的2–6倍.利用長期執行的LASCO拍攝的日冕影象,美國國家航空航天局(NationalAeronauticsandSpaceAdministration,NASA)透過手工記錄的方法建立CoordinatedDataAnalysisWorkshops(CDAW)[1]CME事件庫,但是手動對每個事件標註過於繁瑣且存在個人的主觀偏差.
隨著自動檢測技術的迅速發展,湧現了一系列自動檢測識別CME的方法[2].Robbrecht等[3]基於霍夫變換首次提出一種自動檢測方法ComputerAidedCMETrackingcatalog(CACTus).Olmedo等[4]基於區域增長演算法提出SolarEruptiveEventDetectionSystem(SEEDS).除了以上兩種基於灰度特徵的識別方法,Boursier等[5]提出AutomaticRecognitionofTransientEventsandMarseilleInventoryfromSynopticmaps(ARTEMIS).Goussies等[6]提出了一種基於紋理特徵灰度共生矩陣的非引數監督的CME分割方法.Colaninno等[7]提出了一種基於光流法的CME檢測和跟蹤演算法.Liu等[8]使用支援向量機(SupportVectorMachine,SVM)計算CME到達時間估計Qiang等[9]提出了一種基於自適應背景學習技術檢測CME方法.Zhang等[10]提出了極限學習機(ExtremeLearningMachine,ELM)基於影象亮度和紋理特徵檢測CME,並結合時空連續性排除誤判區域.
以上所述自動檢測方法多為基於灰度特徵、紋理特徵、光流法、傳統的機器學習.由於CME具有多種特徵,這些方法主要基於人為選擇的特徵或利用設定簡單的閾值進行處理,並不能達到很好的檢測效果.而深度學習具有強大的特徵提取功能,自動學習得到有效特徵.Wang等[11]基於卷積神經網路(ConvolutionalNeuraNetwork,CNN)提出了CMEAutomaticdetectionandtrackingwithMachinELearning(CAMEL)自動識別跟蹤CME方法.
隨著大資料和深度學習的發展,CNN在影象分類及計算機視覺領域被廣為使用.通常,CNN使用堆疊的卷積核來逐層提取特徵,每個卷積核僅專注一種特徵.它們在整個影象中共享權重.與全連線的神經網路相比,CNN提高了特徵提取效率,大大減少了計算量,並且可以有效地處理矩陣資料.在太陽活動的分析和研究中,深度學習演算法也引起了天文學家的關注並得到應用[12].Hernandez[13]將卷積神經網路應用於太陽耀斑預測,Huang等[14]採用深度CNN構建太陽耀斑預報模型,Szenicer等[15]使用CNN網路得到極紫外窄帶影象到光譜輻照度測量的對映.Armstrong等[16]基於卷積神經網路的方法,提取SolarOpticalTelescope影象特徵分類為暗條、日珥、耀斑帶、黑子和寧靜太陽.Ahmadzadeh等[17]基於深層網路的方法分割暗條.Wang等[18]使用深度學習框架建立CME到達地球時間的預測模型.
本文采用深層VisualGeometryGroup(VGG)網路,利用LASCOC2的白光日冕儀觀測,對日冕儀影象按照有無觀測到CME進行分類.含有CME的影象標籤為1,反之則標籤為0.此外,基於VGG分類出來的標籤,我們結合了時間序列特性,消除了誤判區域.根據分類結果,我們對CME影象序列進行了時間屬性統計分析,並與CDAW人工事件庫進行了比較.未來先進天基太陽天文臺(AdvancedSpace-basedSolarObservatoryASO-S)[19]衛星的萊曼阿爾法太陽望遠鏡(TheLyman-alphaSolarTelescope,LST)有效載荷上搭載有日冕儀(SolarCoronaImager,SCI)[20,21,22,23].我們將對該儀器產生的日冕影象進行有無CME的分類,標籤為1的影象將推送給國內的各空間天氣預報中心對CME進行預警[24].
2、日冕儀影象分類的深度學習模型方法
本文選取LASCOC2日冕儀6個月的觀測資料,其中2011年1月的影象作為訓練集2011年2月半個月的影象作為測試集,2012年和2014年兩年對應的2月和3月共4個月的影象用於研究分類結果與CDAW比較以及探尋和太陽黑子活動較大較小月份的關係.
2.1 資料預處理
利用SolarSoftware(SSW)中的程式,我們對日冕儀資料進行預處理.使用lascoreadfits.pro讀取0.5級LASCOC2的fits檔案,然後使用reducelevel1.pro將其處理為leve1資料.該處理包括對暗電流、平場、雜散光、畸變、漸暈、輻射定標、時間和位置校正的校準.經過處理後,太陽北已經旋轉到影象北.作為預處理步驟,首先將所有1024×1024畫素的LASCOC2輸入影象降取樣為512×512畫素.然後,所有降取樣的影象都將透過噪聲濾波器,以抑制某些尖銳的噪聲特徵.本文采用了大小為3×3的滑動視窗歸一化塊濾波器.歸一化塊濾波器是一種基本的線性影象濾波器,輸出畫素值是核視窗內畫素值的均值.然後,使用以下公式計算出差分影象:
其中,pt表示當前執行差分影象,nt表示當前影象,nt-1表示上一張影象.
2.2 構建資料集
機器學習主要分有監督學習和無監督學習.有監督學習是指在已知輸入及其對應輸出的情況下,透過訓練這些資料,來發現它們之間的對映關係.無監督學習僅具有輸入資料,而沒有對應的輸出.它需要依靠這些已知資料的特徵統計找到其固有關聯.本文使用有監督學習來解決日冕儀影象的分類問題,檢測影象中是否有CME發生.對預處理完的資料進行標籤分類,從CDAW事件庫中獲取標籤,但是從實際的圖中,我們發現有些圖含有CME結構,而CDAW沒有記錄.因此,在CDAW的基礎上,我們需要再進行人工分類,將2011年1月和2月的資料二次分類.該資料作為本文的訓練集和測試集.
2.3 分類模型
目前,在計算機視覺領域中的深度學習模型為CNN,常用於分類的CNN經典模型有VGG、AlexNet、LeNet[25,26,27],CNN利用影象的空間相關性提取影象的輪廓資訊,提高了網路的學習能力.本文日冕儀影象分類方法採用穩定且高效能的VGG模型.
圖1為VGG16模型結構.首先,本文將預處理完的影象降取樣為224×224畫素作為輸入影象,由於影象為灰度影象,為滿足VGG3通道需求,本文將灰度影象進行復制,分別輸入模型中R、G、B3通道中,將一幅影象表示為224×224×3的矩陣.
圖1基於VGG16的影象分類模型
VGG透過多次堆疊3×3的卷積核和2×2的最大池化層,來構建深層卷積神經網路VGG16有13個卷積層和3個全連線層,其中13個卷積層分別在第2、4、7、10和13層被池化層分割,最大池化層起降維操作、保留最大數值、提高計算速度,同時提高所提取特徵的穩健性.在執行完具有卷積層和池化層的5個迭代過程後,原始的224×224×3特徵圖已縮減為7×7×512.然後執行3個全連線層的操作,7×7×512特徵圖經過第1次全連線操作後的輸出單元為4096,為了減輕和防止過擬合,我們在訓練過程中使用dropout函式先隨機扔掉一部分神經元,再進行第2次全連線操作,該全連線層的輸出也為4096.由於本文為二分類,所以將第3個全連線層的輸出改為2個輸出單元.它們代表了CME發生和未發生的機率,再使用softmax函式進行歸一化計算,求得影象是否有CME結構.
每個卷積層都用3×3的卷積核進行卷積,控制滑動步長,從左到右,從上到下滑動公式可表示為如下:
其中,表示第l層第j個特徵圖,N表示第l-1層特徵圖的數量,表示第l-1層第i個特徵圖,表示第l層第i個特徵圖的卷積核,表示第l層第j個特徵圖的偏差項,f(x)表示非線性啟用函式,max函式表示返回給定引數的最大值.卷積操作之後進入啟用層特徵圖經過非線性啟用函式如sigmoid函式、符號函式(sign)或修正線性單元(RectifiedLinearUnit,ReLU)處理後得到啟用圖.本文使用ReLU函式.將啟用特徵圖再進行最大池化操作.計算每個特徵圖中區域性感受域的最大值,用最大值表示該領域,領域步幅為2在執行完卷積層和全連線層後,使用softmax函式進行分類,公式表示為:
其中,PCME表示測試影象含有CME的機率,xCME和xNOT-CME都是來自最終輸出層的輸出單元.CNN訓練目的是讓損失函式的值達到最小,交叉熵損失公式表示為:
其中,L表示損失值,N表示訓練影象數量,yi表示第i張影象的真實標籤值,ai表示第i張影象softmax求得的預測標籤值.最後我們選擇自適應學習率的Adam最佳化器,Adam帶有動量項的RMSprop,利用梯度的一階矩估計和2階矩估計動態調整每個引數的學習率.
2.4 劃分CME影象序列
我們使用訓練得到的模型,對2012年和2014年兩年對應的2月和3月的影象進行預測,最終得到了預測標籤.如果將連續都是標籤為1的影象歸為一個CME影象序列,有些影象序列是不完整的.因此,結合時空的連續性,需要重新制定規則來分割CME影象序列.首先,允許存在間隔一張圖標籤為0,但不能連續兩張圖標籤為0.按照第1個規則,我們可以得到每個初步劃分的影象序列.接著,對於影象序列的總時間和張數較少的進行進一步操作:丟棄還是保留這個影象序列.如果影象序列的總時間小於0.8h,並且圖片數少於4張,我們丟棄該影象序列.反之則保留該影象序列.最後,對這部分保留下來的影象序列再進行進一步操作:合併到前一個影象序列、合併到後一個影象序列或保留不進行合併.我們分別計算與前後兩個影象序列的時間差,透過設定時間閾值1h來解決.如果與前後影象序列都超過1h,則不進行合併.
3、實驗結果與分析
本文在LASCOC2資料集上進行日冕儀影象分類實驗,使用Pytorch1.2.0框架和Python3.7語言實現,VGG模型在單塊QuadroP5000的GPU上訓練完成.本文選取了2011年1月和2月的影象做訓練集和驗證集,資料集共有4483張影象,其中包括3126張訓練影象和1357張驗證影象.對2011年1月和2月的日冕影象進行降噪等預處理後,輸入到構建好的網路中進行訓練.訓練階段超引數設定:初始學習率為1×10-4,正則化係數為1×10-8,損失函式為CrossEntropyLoss,最佳化器選擇自適應學習率的Adam,模型透過隨機引數初始化開始訓練.訓練完模型後,進入測試階段,本文選取了2012年和2014年兩年對應的2月和3月的日冕儀影象,共12236張影象,進行CME影象序列分類測試.
3.1 模型分析
圖2可看出,經過20輪訓練次數(Epoch),測試集損失(Loss)趨於穩定,測試集在模型上的最高準確率為92.5%.表1為本文使用的VGG模型和Wang等人的LeNet模型[11]得到的分類模型評估比較.計算了準確率(Accuracy)、召回率(Recall)、精準率(Precision)、被模型預測為正的正樣本(TruePositives)、被模型預測為負的負樣本(FalsePositives).可看出本文VGG的分類準確率達到92.5%,高於LeNet模型的86.2%.
3.2 影象分類結果分析
本文總共統計到230個CME影象序列.從圖3中可看出,多數CME影象序列持續時間在2h左右,少數CME影象序列超過5h.我們從圖中可看到時間持續較長的影象序列,最長達到104h.這是因為按標籤以及結合時空連續性產生的CME影象序列中,有些CME事件是連續發生的.部分CME影象序列包含了多個CME事件,進而造成CME影象序列總體持續時間較長.
圖2VGG16模型測試集準確率和損失隨訓練次數的變化
表1與LeNet網路對比結果
圖32012年2、3月和2014年2、3月4個月資料的每個CME影象序列持續時間統計圖.
根據圖4右圖太陽黑子活動年份,本文選取了2012年太陽黑子數較少的2月和3月,2014年太陽黑子數較多的2月和3月,使用箱型圖進行統計分析.箱型圖能提供有關資料位置和分散情況的關鍵資訊,尤其在比較不同的母體資料時更可表現其差異.離群點分佈在箱型圖外側,表現為有些影象序列包含了多個CME事件,導致此類影象序列的總時間很長.另一方面,能夠體現這類影象序列CME活動較劇烈.
圖4左圖為2012年2月至3月(粉色)和2014年2月至3月(藍色)的各CME影象序列的時間統計箱型圖,每個箱型包含5條線,從上至下:上邊緣、上四分位數、中位數、下四分位數、下邊緣,菱形資料點為離群點;右圖為太陽活動黑子數的每天(黃線)、每月(藍線)、每月平滑(紅線)的統計曲線圖,StandardCurve(SC)預測(紅點):僅基於黑子數序列,CombinedMethod(CM)預測(紅破折線):結合黑子數序列和aa地磁指數.
圖4左圖可看出,2012年和2014年這兩年對應的兩個月資料,2014年上四分位數與下四分位數之差較大,而2012年的較小.對比太陽活動年月,2014年2月太陽活動水平高,2012年2月太陽活動水平較低或許與CME活動程度相關.
統計分析本文分類方法篩選出來的每個CME影象序列,與CDAW事件庫比較發現,在起始時刻上,與CDAW記錄的CME事件基本相差在24min內.圖5根據我們的標籤結合時空規則,將圖5的第2張圖至倒數第2張圖結束歸為一個CME影象序列,但CDAW未記錄該時段CME事件,表明本文的模型對CME結構較弱的事件具有較高的靈敏度.其中第2行第3張圖根據2.4節中定義的時空連續性規則,這個單張的標籤為0的影象仍屬於該CME影象序列.
暈狀(halo)CME是和災害性空間天氣最密切相關的一類CME.圖6是2014年2月19日一個暈狀CME事件的影象分類結果.我們發現,該CME影象序列所有圖片的分類標籤全部被成功地標註為1.因而,對於這類較強的CME事件,本文的分類方法具有很高的分類準確率.
圖5CDAW上未標註的CME影象序列舉例,每張圖上方的標籤1表示該影象中含有CME結構,反之標籤0表示不含有CME結構.每張圖下方的T代表時間,時間標準為UniversalTimeCoordinated(UTC).
圖6暈狀CME影象序列舉例
圖7展示一個持續時間較長且含有多個CME事件的CME影象序列.根據分類標籤和時空連續性規則,圖7的第2張至最後一張歸為一個影象序列.從圖中可以發現,該影象序列至少包含了兩個以上的窄型CME事件.對於此類CME影象序列,我們目前還不能將各個CME事件區分開來.CME事件的分離依賴於我們後續的步驟,也就是識別追蹤過程[11].
圖7分類的CME影象序列至少含有兩個CME事件舉例
4、總結與展望
本文選取了部分LASCOC2日冕影象並做預處理,從CDAW事件庫中獲取標籤,但是發現有些圖含有CME結構而CDAW沒有記錄.因此,在CDAW的基礎上,我們再次進行了人工分類.本文使用了VGG16卷積神經網路模型,同時結合時空連續性規則,能夠自動有效分類出各種CME影象序列,甚至檢測出較弱的CME結構.測試集影象分類準確率達到92.5%,優於Wang等[11]檢測CME使用的LeNet模型結果.對於CME活動較劇烈的時間段,分類出的一個CME影象序列可能包含有至少兩個CME事件.與CDAW事件庫比較,本文分類出的CME影象序列包含了絕大部分CDAW標註的CME事件.在CME開始發生時刻上,本文與CDAW標註的時間基本相差在24min內.後續我們將統計分析更多的LASCOC2影象資料,並對CME進行識別和檢測跟蹤來提取各個CME的主要引數並建立資料庫.未來本文的方法將應用到ASO-S衛星上,對SCI產生的日冕影象進行有無CME結構的影象分類,建立CME標籤庫,推送給合作的空間天氣預報中心,對CME到達地球的時間進行預報.
