基于深度学习的极光亚暴时-空自动检测

杨秋菊，任杰，向晗

陕西师范大学物理与信息技术学院,西安 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提供.在此感谢上述机构和专家对本文工作的帮助.