基于机器学习的日冕仪图像分类方法研究∗

单家辉 封 莉 袁汉青 张 岩‡ 钟 忺 甘为群黎辉 黄宇

(1中国科学院紫金山天文台南京210023)

(2中国科学院暗物质和空间天文重点实验室南京210023)

(3武汉理工大学计算机科学与技术学院武汉430070)

(4中国科学技术大学天文与空间科学学院合肥230026)

(5南京大学计算机软件新技术国家重点实验室南京210023)

1 引言

日冕物质抛射(Coronal Mass Ejection,CME)是太阳大气中最剧烈、尺度最大的活动现象,表现为在短时间内日冕结构发生明显的变化,并伴有1011–1013kg携带磁场的等离子体抛射.当日冕物质抛射的方向朝着地球时,可能会与地球磁层发生相互作用,引起近地空间的地磁暴、极光等现象,会对通讯系统和电力系统等产生干扰,严重时会造成巨大的经济损失.因此,CME到达地球的实时预报对空间天气环境的监测十分重要.

CME的自动标注和检测是实现CME预报的重要前提.太阳和日球层天文台(Solar and Heliospheric Observatory,SOHO)搭载的大角度光谱日冕仪(Large Angle and Spectrometric Coronagraph Experiment,LASCO)能够观测太阳日冕活动.LASCO由3台视场不同的日冕仪构成,其中LASCO C2视场的范围大约是太阳直径的2–6倍.利用长期运行的LASCO拍摄的日冕图像,美国国家航空航天局(National Aeronautics and Space Administration,NASA)通过手工记录的方法建立Coordinated Data Analysis Workshops(CDAW)[1]CME事件库,但是手动对每个事件标注过于繁琐且存在个人的主观偏差.随着自动检测技术的迅速发展,涌现了一系列自动检测识别CME的方法[2].Robbrecht等[3]基于霍夫变换首次提出一种自动检测方法Computer Aided CME Tracking catalog(CACTus).Olmedo等[4]基于区域增长算法提出Solar Eruptive Event Detection System(SEEDS).除了以上两种基于灰度特征的识别方法,Boursier等[5]提出Automatic Recognition of Transient Events and Marseille Inventory from Synoptic maps(ARTEMIS).Goussies等[6]提出了一种基于纹理特征灰度共生矩阵的非参数监督的CME分割方法.Colaninno等[7]提出了一种基于光流法的CME检测和跟踪算法.Liu等[8]使用支持向量机(Support Vector Machine,SVM)计算CME到达时间估计.Qiang等[9]提出了一种基于自适应背景学习技术检测CME方法.Zhang等[10]提出了极限学习机(Extreme Learning Machine,ELM)基于图像亮度和纹理特征检测CME,并结合时空连续性排除误判区域.以上所述自动检测方法多为基于灰度特征、纹理特征、光流法、传统的机器学习.由于CME具有多种特征,这些方法主要基于人为选择的特征或利用设定简单的阈值进行处理,并不能达到很好的检测效果.而深度学习具有强大的特征提取功能,自动学习得到有效特征.Wang等[11]基于卷积神经网络(Convolutional Neural Network,CNN)提出了CME Automatic detection and tracking with MachinE Learning(CAMEL)自动识别跟踪CME方法.

随着大数据和深度学习的发展,CNN在图像分类及计算机视觉领域被广为使用.通常,CNN使用堆叠的卷积核来逐层提取特征,每个卷积核仅专注一种特征.它们在整个图像中共享权重.与全连接的神经网络相比,CNN提高了特征提取效率,大大减少了计算量,并且可以有效地处理矩阵数据.在太阳活动的分析和研究中,深度学习算法也引起了天文学家的关注并得到应用[12].Hernandez[13]将卷积神经网络应用于太阳耀斑预测,Huang等[14]采用深度CNN构建太阳耀斑预报模型,Szenicer等[15]使用CNN网络得到极紫外窄带图像到光谱辐照度测量的映射.Armstrong等[16]基于卷积神经网络的方法,提取Solar Optical Telescope图像特征分类为暗条、日珥、耀斑带、黑子和宁静太阳.Ahmadzadeh等[17]基于深层网络的方法分割暗条.Wang等[18]使用深度学习框架建立CME到达地球时间的预测模型.

本文采用深层Visual Geometry Group(VGG)网络,利用LASCO C2的白光日冕仪观测,对日冕仪图像按照有无观测到CME进行分类.含有CME的图像标签为1,反之则标签为0.此外,基于VGG分类出来的标签,我们结合了时间序列特性,消除了误判区域.根据分类结果,我们对CME图像序列进行了时间属性统计分析,并与CDAW人工事件库进行了比较.未来先进天基太阳天文台(Advanced Space-based Solar Observatory,ASO-S)[19]卫星的莱曼阿尔法太阳望远镜(The Lyman-alpha Solar Telescope,LST)有效载荷上搭载有日冕仪(Solar Corona Imager,SCI)[20–23].我们将对该仪器产生的日冕图像进行有无CME的分类,标签为1的图像将推送给国内的各空间天气预报中心,对CME进行预警[24].

2 日冕仪图像分类的深度学习模型方法

本文选取LASCO C2日冕仪6个月的观测数据,其中2011年1月的图像作为训练集,2011年2月半个月的图像作为测试集,2012年和2014年两年对应的2月和3月共4个月的图像用于研究分类结果与CDAW比较以及探寻和太阳黑子活动较大较小月份的关系.

2.1 数据预处理

利用Solar Software(SSW)中的程序,我们对日冕仪数据进行预处理.使用lasco_readf its.pro读取0.5级LASCO C2的f its文件,然后使用reduce_level_1.pro将其处理为level 1数据.该处理包括对暗电流、平场、杂散光、畸变、渐晕、辐射定标、时间和位置校正的校准.经过处理后,太阳北已经旋转到图像北.作为预处理步骤,首先将所有1024×1024像素的LASCO C2输入图像降采样为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–27],CNN利用图像的空间相关性提取图像的轮廓信息,提高了网络的学习能力.本文日冕仪图像分类方法采用稳定且高性能的VGG模型.

图1为VGG 16模型结构.首先,本文将预处理完的图像降采样为224×224像素作为输入图像,由于图像为灰度图像,为满足VGG 3通道需求,本文将灰度图像进行复制,分别输入模型中R、G、B 3通道中,将一幅图像表示为224×224×3的矩阵.

图1 基于VGG 16的图像分类模型Fig.1 Image classif ication model based on VGG 16

VGG通过多次堆叠3×3的卷积核和2×2的最大池化层,来构建深层卷积神经网络.VGG 16有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的卷积核进行卷积,控制滑动步长,从左到右,从上到下滑动,公式可表示为如下:

其中,PCME表示测试图像含有CME的概率,xCME和xNOT-CME都是来自最终输出层的输出单元.CNN训练目的是让损失函数的值达到最小,交叉熵损失公式表示为:

其中,L表示损失值,N表示训练图像数量,yi表示第i张图像的真实标签值,ai表示第i张图像softmax求得的预测标签值.最后我们选择自适应学习率的Adam优化器,Adam带有动量项的RMSprop,利用梯度的一阶矩估计和2阶矩估计动态调整每个参数的学习率.

2.4 划分CM E图像序列

我们使用训练得到的模型,对2012年和2014年两年对应的2月和3月的图像进行预测,最终得到了预测标签.如果将连续都是标签为1的图像归为一个CME图像序列,有些图像序列是不完整的.因此,结合时空的连续性,需要重新制定规则来分割CME图像序列.首先,允许存在间隔一张图标签为0,但不能连续两张图标签为0.按照第1个规则,我们可以得到每个初步划分的图像序列.接着,对于图像序列的总时间和张数较少的进行进一步操作:丢弃还是保留这个图像序列.如果图像序列的总时间小于0.8 h,并且图片数少于4张,我们丢弃该图像序列.反之则保留该图像序列.最后,对这部分保留下来的图像序列再进行进一步操作:合并到前一个图像序列、合并到后一个图像序列或保留不进行合并.我们分别计算与前后两个图像序列的时间差,通过设定时间阈值1 h来解决.如果与前后图像序列都超过1 h,则不进行合并.

3 实验结果与分析

本文在LASCO C2数据集上进行日冕仪图像分类实验,使用Pytorch 1.2.0框架和Python 3.7语言实现,VGG模型在单块Quadro P5000的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)、被模型预测为正的正样本(True Positives)、被模型预测为负的负样本(False Positives).可看出本文VGG的分类准确率达到92.5%,高于LeNet模型的86.2%.

3.2 图像分类结果分析

本文总共统计到230个CME图像序列.从图3中可看出,多数CME图像序列持续时间在2 h左右,少数CME图像序列超过5 h.我们从图中可看到时间持续较长的图像序列,最长达到104 h.这是因为按标签以及结合时空连续性产生的CME图像序列中,有些CME事件是连续发生的.部分CME图像序列包含了多个CME事件,进而造成CME图像序列总体持续时间较长.

图2 VGG 16模型测试集准确率和损失随训练次数的变化Fig.2 VGG 16 model test set accuracy and loss change with ep och

表1 与LeNet网络对比结果Table 1 Comp arison with the r esults of LeNet

图3 2012年2、3月和2014年2、3月4个月数据的每个CME图像序列持续时间统计图.Fig.3 Statistics of the duration of each CME image sequence in February and March 2012,and February and M arch 2014.

根据图4右图太阳黑子活动年份,本文选取了2012年太阳黑子数较少的2月和3月,2014年太阳黑子数较多的2月和3月,使用箱型图进行统计分析.箱型图能提供有关数据位置和分散情况的关键信息,尤其在比较不同的母体数据时更可表现其差异.离群点分布在箱型图外侧,表现为有些图像序列包含了多个CME事件,导致此类图像序列的总时间很长.另一方面,能够体现这类图像序列CME活动较剧烈.

图4 左图为2012年2月至3月(粉色)和2014年2月至3月(蓝色)的各CME图像序列的时间统计箱型图,每个箱型包含5条线,从上至下:上边缘、上四分位数、中位数、下四分位数、下边缘,菱形数据点为离群点;右图为太阳活动黑子数的每天(黄线)、每月(蓝线)、每月平滑(红线)的统计曲线图,Standard Curve(SC)预测(红点):仅基于黑子数序列,Combined Method(CM)预测(红破折线):结合黑子数序列和aa地磁指数.图源http://www.sidc.be/images/wolfjmms.png.Fig.4 The left panel is a box diagram of time statistics of each CME image sequence from February to M arch 2012(pink)and February to March 2014(blue).Each box contains f ive lines,which from top to bottom indicate respectively the upp er edge,upper quantile,median,lower quartile,and lower edge.And diamond data p oints are outlier.The right panel is the statistical curve of the number of sunsp ots p er day(yellow line),p er month(blue line),and p er month after smoothing(red line).SC Predictions(red dots):only based on the sunsp ot number series.CM Predictions(red dashes):combining the sunspot number series with the aa geomagnetic index.Image source http://www.sidc.be/images/wolfjmms.png.

图4左图可看出,2012年和2014年这两年对应的两个月数据,2014年上四分位数与下四分位数之差较大,而2012年的较小.对比太阳活动年月,2014年2月太阳活动水平高,2012年2月太阳活动水平较低或许与CME活动程度相关.

统计分析本文分类方法筛选出来的每个CME图像序列,与CDAW事件库比较发现,在起始时刻上,与CDAW记录的CME事件基本相差在24 min内.图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事件,本文的分类方法具有很高的分类准确率.

图5 CDAW上未标注的CME图像序列举例,每张图上方的标签1表示该图像中含有CME结构,反之标签0表示不含有CME结构.每张图下方的T代表时间,时间标准为Universal Time Coordinated(UTC).Fig.5 An example of the CM E image sequence which is not catalogued by CDAW.In each panel,label 1 indicates an image with CME structures,while label 0 without.The T at the bottom of each panel represents time,and the time standard is UTC.

图6 晕状CME图像序列举例Fig.6 An example of a halo CME image sequence

图7 展示一个持续时间较长且含有多个CME事件的CME图像序列.根据分类标签和时空连续性规则,图7的第2张至最后一张归为一个图像序列.从图中可以发现,该图像序列至少包含了两个以上的窄型CME事件.对于此类CME图像序列,我们目前还不能将各个CME事件区分开来.CME事件的分离依赖于我们后续的步骤,也就是识别追踪过程[11].

图7 分类的CME图像序列至少含有两个CME事件举例Fig.7 An example of a CME image sequence containing at least two CME events

4 总结与展望

本文选取了部分LASCO C2日冕图像并做预处理,从CDAW事件库中获取标签,但是发现有些图含有CME结构而CDAW没有记录.因此,在CDAW的基础上,我们再次进行了人工分类.本文使用了VGG 16卷积神经网络模型,同时结合时空连续性规则,能够自动有效分类出各种CME图像序列,甚至检测出较弱的CME结构.测试集图像分类准确率达到92.5%,优于Wang等[11]检测CME使用的LeNet模型结果.对于CME活动较剧烈的时间段,分类出的一个CME图像序列可能包含有至少两个CME事件.与CDAW事件库比较,本文分类出的CME图像序列包含了绝大部分CDAW标注的CME事件.在CME开始发生时刻上,本文与CDAW标注的时间基本相差在24 min内.后续我们将统计分析更多的LASCO C2图像数据,并对CME进行识别和检测跟踪来提取各个CME的主要参数并建立数据库.未来本文的方法将应用到ASO-S卫星上,对SCI产生的日冕图像进行有无CME结构的图像分类,建立CME标签库,推送给合作的空间天气预报中心,对CME到达地球的时间进行预报.