基于深度學習的極光亞暴時-空自動檢測

楊秋菊，任杰，向晗

陜西師范大學物理與信息技術學院,西安 710061

0 引言

亞暴是一種發生在地球夜半球太空區域的巨大能量瞬間釋放現象(Akasofu,1964),亞暴發生過程伴隨著極光形態和亮度的劇烈變化,研究人員將這種極光活動稱為極光亞暴.極光亞暴的發生伴隨著由電離層電流引起的強烈磁干擾(Akasofu,2017),會對人類活動造成很大的影響,例如引起地球同步軌道上的衛星充電(李柳元等,2006)、高緯度地區無線電通訊中斷(Wen et al.,2012;Chernyshov et al.,2020)和供電系統故障(Boteler,2019)等.目前,極光亞暴的產生機制仍存在爭議(Schindler,1974;Angelopoulos et al.,2008;Maimaiti et al.,2019),行星際磁場由北轉向南是否是亞暴發生的原因仍有待確認.從海量的觀測數據中自動檢測亞暴事件并確定其發生位置有助于統計分析亞暴的產生機制(Ieda et al.,2019).

目前,檢測極光亞暴事件的途徑可以分為三種.第一種途徑是從海量極光圖像中人工挑選亞暴事件,如Liou(2010)基于Polar衛星的全域極光圖像手動標記極光亞暴起始時刻及其發生位置(磁地方時-地磁緯度),Ieda等(2019)結合全域極光圖像和全天空極光圖像確定亞暴起始時刻,但人工標注費時費力且在圖像數據量非常大時容易出現漏標、誤標等問題(He et al.,2014).第二種途徑是通過各類物理參數自動檢測亞暴事件,如Sutcliffe(1997)基于Pi2脈沖利用神經網絡的方法自動檢測亞暴事件,Maimaiti等(2019)結合太陽風速度、質子數密度和行星際磁場基于深度學習的方法自動檢測亞暴.雖然該途徑克服了人工挑選的局限性,但亞暴事件真實發生時間和所測量的物理指標之間存在時延(Meng and Liou,2004;Ieda et al.,2019),導致檢測結果存在一定的誤差.為了解決以上兩種問題,人們提出了第三種途徑,基于衛星紫外成像儀(Ultraviolet Imager,UVI)圖像利用機器學習方法自動檢測亞暴事件.如楊秋菊等(2013)基于模糊C均值聚類通過分析亮斑的變化情況確定亞暴的起始時刻,該方法雖然實現簡單并取得了較高的召回率,但準確率有待提升且需要人工設置多個閾值得到檢測結果;結合亞暴事件發生區域的地磁先驗信息,利用SCSLD(Shape-Constrained Sparse and Low-Rank Decomposition)算法通過粗檢測、序列運動分析、細檢測實現亞暴事件檢測,該方法在多年的極光觀測上獲得了很好的檢測效果(Yang et al.,2016),但實現過程中需要人工多次參與分析算法結果,操作復雜而且運行效率較低,實際應用難度大.

基于深度學習的視頻時序行為檢測研究(Zhao et al.,2017;Lin et al.,2018)為亞暴自動檢測提供了新的選擇,但與現有時序行為檢測不同的是亞暴事件檢測重點關注極光點亮的時刻(onset),準確檢測亞暴事件需要構建密集的特征序列;此外,衛星觀測數據時間分辨率較低(0.0056幀/s VS常規視頻24幀/s),極光亞暴持續的視頻幀長度較短,模型需要包含輸入圖像序列的全部時序信息.基于此,本文利用深度學習技術提出了一種端到端的亞暴檢測網絡(Substorm Onset Detection Network,SODN),該方法在取得亞暴檢測高準確率的同時,極大提高了亞暴檢測速度.SODN由兩部分組成:(1)亞暴特征提取模塊;(2)亞暴起始檢測模塊.亞暴特征提取模塊基于雙流網絡(Simonyan and Zisserman,2014)構建,包含兩個并行的卷積神經網絡:空間網絡和時序網絡.其中,空間網絡用于提取亞暴空間形態特征,時序網絡用于提取亞暴時序運動特征,本文將上述兩種特征級聯獲得亞暴特征表示.亞暴起始檢測模塊利用三個一維時序卷積(Lea et al.,2016)融合上述兩種特征,輸出亞暴起始的概率序列.與第三種途徑不同的是本文方法利用雙流網絡自動提取亞暴時-空特征,不需要手動設計亞暴特征表示.此外,利用三個一維時序卷積層構建亞暴起始檢測模塊獲取亞暴起始的邊界信息,無需繁瑣的篩選分析過程直接得出檢測結果,極大提高了亞暴檢測效率.

本文的結構安排如下:第1節簡要介紹極光亞暴事件特點和實驗數據來源,第2節詳細介紹SODN框架,實驗結果和分析在第3節給出,第4節為結論和討論.

1 數據集與數據預處理

1.1 UVI圖像

本文實驗數據由Polar衛星攜帶的UVI拍攝(Torr et al.,1995).該衛星于1996年2月24日進入高橢圓形軌道,遠地點高度為9個地球半徑,近地點高度1.5個地球半徑,軌道傾角86°,周期約17.5 h.Polar衛星于2008年4月完成相關計劃,在此期間收集了大量的UVI極光圖像,首次詳細觀察了地磁亞暴發生期間極光的演變過程.Polar衛星的軌道平面在觀測期間緩慢向南移動,本文選用Polar衛星位于北半球期間1996年12月,1997年1月、2月、12月和1998年12月采集的UVI圖像作為數據集,其中共包含371個亞暴事件;從每月中隨機選取5天作為測試數據集,共包含亞暴事件66個;其余為訓練數據集,共包含亞暴事件305個.

UVI采集了多個波段的極光圖像(Torr et al.,1995),借鑒文獻(Kullen and Karlsson,2004)的處理方式本文只使用LBHL(Lyman-Birge-Hopfield long)濾波通道下的UVI圖像.該波段下UVI圖像的空間分辨率為30 km,圖像覆蓋范圍約4104000 km2,時間分辨率為3 min,大小為228×200;UVI圖像信息除了像素亮度外,還包含了每個像素對應的MLT、MLAT等信息(楊秋菊等,2013).本文采用Liou(2010)提供的人工標注亞暴事件列表,將該標注作為挑選訓練數據集和衡量亞暴檢測結果的基準.

1.2 數據預處理

Polar衛星采集UVI圖像時位置在不斷變化,UVI圖像的同一位置在不同圖像中對應著不同的磁地方時和地磁緯度(楊秋菊等,2013),不利于從圖像序列中統計亞暴發生時亮斑的位置變化.本文使用和我們前期工作(楊秋菊等,2013)相同的方法將UVI圖像轉換至磁地方時-地磁緯度(MLT-MLAT)坐標下,轉換后圖像中的相同像素位置對應相同的磁地方時和地磁緯度,具體步驟為:

(1)根據亞暴發生的時間和位置,重點關注夜側和高緯區域,MLT-MLAT網格圖像限定為磁地方時18—24 MLT及00—06 MLT,地磁緯度50°—90°MLAT;

(2)設定網格大小,轉換后的網格圖像大小為200×120((90°—50°MLAT)/0.2°=200行,12 h/0.1 h=120列);

(3)對網格中每個坐標點,找到其在UVI圖像中的位置,計算UVI圖像中落入每個網格的像素的均值,該均值作為每個網格點的像素強度.

由MLT-MLAT網格圖像的定義知,各網格圖像相同位置對應相同的MLT和MLAT值,極光亞暴過程中亮斑區域的形態變化在MLT-MLAT網格圖像中主要表現為水平方向和豎直方向的形態變化(如圖1所示),極光形態隨時間的演變規律可以直接從MLT-MLAT圖像序列中計算.

圖1為一組轉換后的MLT-MLAT網格圖像示例，該圖像展示了一個典型的極光亞暴事件發生過程,根據該過程中極光的形態特征可以將亞暴事件分為三個階段(Rostoker et al.,1980;Liang et al.,2018):

圖1 1996年12月04日21∶51∶41—22∶28∶10UT采集的UVI圖像轉為MLT-MLAT網格圖像的結果，MLT范圍為18—06，MLAT范圍為50°—90°.圖中粗體標注的時間為亞暴膨脹階段，其中22∶00∶34UT為亞暴事件的起始時刻，斜體標注的時間為亞暴消退階段Fig.1 The MLT-MLAT grid images of the UVI images captured on December 04,1996 between 21∶51∶41—22∶28∶10UT.The range of MLT is 18—06,and the range of MLAT is 50°—90°.The time marked in bold is the substorm expansion phase,where 22∶00∶34UT is the onset of the substorm,and the time marked in italics is the substorm recovery phase

(1)增長階段:地面成像儀可以觀測到極光弧向赤道向移動,衛星成像儀觀測到極光卵區域仍保持靜止.

(2)膨脹階段:以極光卵局部區域突然點亮為開始標志,點亮區域迅速向極向和東西向膨脹,同時伴隨著亮度的增強,最后到達一個極向最大邊界.

(3)消退階段:該階段發生在點亮區域膨脹至極向最大位置后,表現為點亮區域亮度逐漸減弱、膨脹區域消退,最后恢復到亞暴之前的水平.

其中,增長階段的持續時間通常為30～60 min,膨脹階段和消退階段共持續1～2 h(Ebihara,2019).膨脹期內,從UVI圖像上觀測亞暴的兩個關鍵特征(Kepko and McPherron,2001;Frey,2004)是:(1)極光卵部分區域極光亮度迅速增強;(2)亮度增強區域向極向和赤道方向膨脹,尤其在磁地方時(MLT)21—03時和地磁緯度(MLAT)65°—75°之間.只有空間形態變化和時序變化同時滿足上述兩個特征才被稱為亞暴,只滿足第一個特征可能是偽暴(Nakamura et al.,1994;Aikio et al.,1999;Kullen and Karlsson,2004)或圖像中的噪聲點.基于此,SODN使用空間網絡和時序網絡分別關注亞暴形態信息和時序運動信息,最后通過時序卷積融合兩種信息獲得檢測結果.

2 SODN模型

本文提出的SODN模型如圖2所示.首先將轉換后的MLT-MLAT圖像序列分為若干個連續等長的片段,每個片段生成一個單元,一個單元包含片段內的單幀圖像和片段內的堆疊光流圖(x方向和y方向).然后利用特征提取模塊分別對每個單元提取亞暴空間特征和時序特征,并級聯兩個特征作為單元級特征向量.最后按時間順序拼接每個單元的特征向量作為亞暴起始檢測模塊的輸入,輸出為整個UVI圖像序列的亞暴起始概率向量,向量中的每個值為其對應單元作為亞暴起始的置信度分數.

2.1 亞暴特征提取模塊

為了提取亞暴時-空特征,本文采用雙流網絡構建亞暴特征提取模塊,該網絡由空間網絡和時序網絡組成.如圖2所示,本文將單幀MLT-MLAT圖像輸入空間網絡提取亞暴空間形態特征,將堆疊的光流圖輸入時序網絡提取亞暴時序運動特征.光流計算采用TV-L1算法(Zach et al.,2007),空間網絡和時序網絡都采用具有多個歸一化層的BN-Inception(Ioffe and Szegedy,2015)作為其骨干網絡.該網絡由多個Inception模塊組成,在Inception網絡的基礎上使用多個小尺度卷積代替尺度較大的卷積降低模型參數量,并通過增加BN層使得模型更快收斂,使該模型在準確性和效率之間實現了較好的平衡(Wang et al.,2016).本文在BN-Inception的最后一個全連接層之前加入一個維度為50的全連接層,用于輸出空間特征向量和時序特征向量.

圖2 SODN模型示意圖Fig.2 The diagram of SODN

2.2 亞暴起始檢測模塊

圖3 亞暴起始時刻檢測模塊輸出的概率序列示意圖.選擇輸出概率值大于0.5的時刻作為候選亞暴起始點,在候選起始點中選擇概率峰值點作為檢測結果Fig.3 The diagram of the probability sequence outputted by the substorm onset detection module.The moments when the output probability is greater than 0.5 are selected as the candidate substorm onset points,and we choose the time corresponding to the maximum probability as the final detection result

一個時序卷積層可以表示為：Conv1d(cf,ck,Act)，其中Conv1d表示一維卷積操作，cf,ck,Act分別表示一維卷積的卷積核數量、卷積核尺寸和卷積層的激活函數.亞暴起始檢測器由三個連續的時序卷積層組成，定義為:

Conv1d(256,3,Relu)→Conv1d(256,3,Relu)

→Conv1d(1,3,Sigmoid),(1)

3 實驗結果及分析

3.1 SODN和SCSLD方法的結果對比

為了充分評估SODN與SCSLD的檢測性能,本文選用準確率(precision)、召回率(recall)、F1值和每秒幀數(Frame Per Second,FPS)來衡量兩者的檢測精度和檢測效率,

(2)

(3)

(4)

其中，TP(true positive)表示檢測的亞暴事件在人工標記的亞暴事件列表中;FP(false positive)表示檢測的亞暴事件未出現在人工標記的亞暴事件列表中;FN(false negative)表示出現在人工標記的亞暴事件列表中的事件未被檢測到.

SODN與SCSLD的對比結果如表1所示.SCSLD方法在檢測亞暴過程中,人工多次參與了結果分析和篩選,所以準確率比SODN高1.83%.但從檢測效率方面來看，SODN處理一幀圖像僅需0.003 s，檢測速度高達393幀/s；更重要的是，將UVI圖像輸入SODN后能直接輸出檢測結果，即SODN實現的是端到端的亞暴自動檢測，在實際應用中非常簡便.而相比之下，SCSLD方法整個檢測流程復雜，需要經過粗檢測、時序運動分析、細檢測三步完成檢測，且在時序運動分析中需要對圖像進行兩次運動特征提取，僅其中一次運動特征提取就需要0.15 s(6.93幀/s)，完成檢測一幀的時間很長，從而使其檢測速度很低.

表1 SODN與SCSLD檢測結果Table 1 The detection results of SODN and SCSLD

3.2 SODN與人工標記的結果對比

SODN與人工標記的結果對比如表2所示，SODN從測試數據集共檢測出64個亞暴事件，其中與人工標記(Liou,2010)一致的56個，多檢亞暴事件8個,漏檢亞暴事件10個.分析SODN檢測結果中與人工標記不一致的案例(多檢和漏檢),我們根據其原因不同將這些案例大致歸為以下三類.

表2 SODN與人工標記檢測結果Table 2 The detection results of SODN and manual approach

(1)本文僅考慮了LBHL波段的觀測數據，丟失了亞暴起始關鍵信息

如圖4所示,人工標注列表中1996年12月11日09∶53∶53—09∶54∶36UT為亞暴事件起始時段,但SODN漏檢了該事件.這是因為人工標注亞暴事件時參考了多個波段的UVI極光圖像,而本文僅考慮了LBHL波段,對于圖4所示的亞暴事件來說,其演化速度較快,LBHL波段未采集上亞暴膨脹期極光亮斑從點亮到膨脹的關鍵信息,導致該事件從我們的數據集上看更像一個偽暴事件.類似的情況在1997年12月03日16∶16∶08—16∶16∶55UT也有出現.

圖4 1996年12月11日09∶53∶36—10∶05∶52UT采集的UVI圖像，MLT范圍為18—06，MLAT范圍為50°—90°Fig.4 The UVI images captured on December 11,1996 between 09∶53∶36—10∶05∶52UT.The range of MLT is 18—06,and the range of MLAT is 50°—90°

(2)人工標注中漏標了部分真實存在的亞暴事件

圖5所示為1996年12月14日17∶46∶48—19∶32∶17UT采集的UVI圖像序列.分析該時段UVI圖像序列,用粗體標注的多幀圖像符合亞暴點亮和膨脹相的特征,是一個典型的亞暴事件,但在人工標注列表中該事件未被標注.該天13∶26∶44UT也存在類似的情況,而這些事件均被SODN檢測為亞暴事件.

圖5 1996年12月14日17∶46∶48—18∶14∶24UT采集的UVI圖像序列，MLT范圍為18—06，MLAT范圍為50°—90°Fig.5 The UVI images captured on December 14,1996 between 17∶46∶48—18∶14∶24UT.The range of MLT is 18—06,and the range of MLAT is 50°—90°

(3)較為復雜的亞暴事件

圖6所示為1998年12月14日14∶28∶39—14∶29∶16UT采集的UVI圖像序列,人工標注該時段為亞暴起始時段.分析該時段UVI圖像序列,14∶29∶36UT對應UVI圖像中白色圓圈標記一塊亮斑,其上方存在兩塊細長條亮斑且亮度相近,隨后該亮斑出現微弱的亮度增強和膨脹現象.上述情況增加了該亞暴事件的檢測難度且訓練數據集中類似的亞暴事件很少,導致了SODN漏檢了該事件.

圖6 1998年12月14日14∶29∶36—14∶31∶26UT采集圖像序列，MLT范圍為18—06，MLAT范圍為50°—90°Fig.6 The UVI images captured on December 14,1998 between 14∶29∶36—14∶31∶26UT.The range of MLT is 18—06,and the range of MLAT is 50°—90°

3.3 空間定位統計

亞暴發生的位置信息對亞暴研究同樣具有重要意義.Liou(2010)利用區域生長算法確定亞暴起始時刻圖像中的亮斑區域,再計算亮斑的質心位置作為亞暴發生的位置(MLT-MLAT坐標)進行統計分析.該方法的缺點是每處理一張圖像都需要采用人機交互的方式確認區域生長的種子點.為了解決這一問題,本文利用Grad-CAM(Selvaraju et al.,2017)方法定位亞暴發生位置.該方法利用網絡輸出對指定的卷積層輸出計算梯度得到一個二維矩陣,矩陣中的每個數值表示該區域對網絡決策分類的重要程度,將該二維矩陣渲染為偽彩圖與輸入圖像疊加生成類激活圖(Class Activation Map，CAM)，從圖中可以直觀看出網絡更為關注的區域.本文選取SODN空間網絡的卷積層輸出結果計算CAM圖,利用該圖自動定位亞暴起始時刻亮斑區域,將定位的亮斑區域作為亞暴發生位置進行統計分析.雖然該定位區域可能無法與整個亮斑完全重合,但該區域是整個亮斑區域的子集,其統計結果仍可揭示亞暴發生位置的分布規律.定位過程如圖7所示,具體步驟如下:

圖7 空間網絡(BNInception V2)中Inception5生成CAM圖.從左到右的每一列依次表示MLT-MLAT圖像、生成的CAM圖、通過CAM圖粗略定位亮斑位置、截取出的圖像亮斑部分和最終定位到的亮斑區域.(a)圖中的UVI圖像發生于1998年12月04日12∶12∶49UT,(b)圖中的UVI圖像發生于1998年12月11日08∶13∶45UTFig.7 CAM map generated by Inception5 module in spatial network (BNInception V2).Each column from left to right represents the MLT-MLAT image,the generated CAM map,the rough spot location map through the CAM map,the extracted bright spot,and the final detected spot locations respectively.The UVI image in (a)was captured at 12∶12∶49UT on December 04,1998,and the UVI image in (b)was captured at 08∶13∶45UT on December 11,1998

(1)利用SODN時序檢測結果確定亞暴事件起始幀;

(2)將起始幀輸入空間網絡,獲得空間網絡最后一個Inception模塊的CAM圖;

(3)通過對CAM圖設置閾值(設置多組候選值篩選后得到)獲得圖中權重較大的區域，該區域大致定位亮斑所在位置;

(4)對定位出的區域繪制其灰度直方圖,由于亮斑區域相對其他區域占比較大且亮度高,在直方圖上呈現為一個峰值，將小于峰值對應的灰度值范圍灰度置零，剩下的非零區域就是亮斑所在位置.

利用上述方法對所有檢測出的亞暴事件統計分析其發生位置.圖8給出了亞暴發生位置在MLT-MLAT坐標下的分布情況.圖中亞暴發生位置集中在20—02 MLT、60°—70°MLAT之間、該結果與以往亞暴研究的物理結論高度一致(Liou,2010).

圖8 亞暴起始時刻亮斑在MLT-MLAT網格圖中的分布位置及發生頻率(%)Fig.8 The distribution and occurrence frequency (%)of the substorm onset locations in MLT-MLAT grid chart

4 結論

針對現有亞暴自動檢測方法的不足,本文提出一種基于深度學習的亞暴事件時-空自動檢測模型SODN.該模型利用兩個并行的卷積神經網絡分別提取亞暴的空間形態特征和時序運動特征,然后將兩種特征級聯生成亞暴特征表示.該特征可以有效表示亞暴從增長相到恢復相整個過程中極光卵在時空上的變化情況,解決了傳統方法需要分步提取亞暴特征的問題,簡化了特征提取流程，降低了模型復雜度,使得該模型在實際應用中擁有良好的部署性和擴展性.同時為了兼顧亞暴檢測的準確性和效率,SODN利用三個緊湊的時序卷積層融合多個時序范圍的亞暴特征,直接輸出亞暴起始的概率序列,克服了傳統自動檢測方法需要人工參與分析、手動設置閾值和檢測效率較低的問題.在多年的Polar衛星觀測數據上SODN實現了端到端的亞暴事件自動檢測,檢測速度高達393幀/s;而且利用該方法對亞暴發生位置的統計結果與已有物理結論非常吻合,證明了SODN可以應用于大規模亞暴時-空檢測,具有良好的實用價值.

亞暴特征提取的優劣會直接影響亞暴檢測結果,后續工作將嘗試采用最新的自監督光流模型(Sun et al.,2018)改進光流效果,進一步提升檢測準確性.此外,數據集的質量直接影響著網絡的性能,本文僅考慮單個波段的數據，時間分辨率較低且部分數據存在標注錯誤,因此SODN在檢測性能上仍有較大的提升空間.后續工作將結合衛星多個波段采集的UVI圖像提升數據集的時間分辨率,并將該模型應用于其他衛星(如IMAGE衛星)觀測數據,進一步驗證模型的泛化性能.

致謝本文UVI圖像數據由美國國家空間科學和技術中心(National Space Science and Technology Center,NSSTC)提供,亞暴事件的人工標記由Liou提供.在此感謝上述機構和專家對本文工作的幫助.