机器学习基础上日冕仪图像分类方法探究

　　摘要：日冕物质抛射(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到达地球的时间进行预报.