梁宇峰，白先勇,2,3，冯 松

(1. 昆明理工大学 云南省计算机技术应用重点实验室/信息工程与自动化学院,云南 昆明 650500；2.中国科学院 太阳活动重点实验室,北京100012；3.中国科学院 国家天文台,北京 100049)

太阳活动是太阳大气中局部区域的各种不同活动现象的总称。当太阳活动增强时会突然释放出巨大的能量，同时抛射出不同能量的粒子对地球空间环境造成很大的影响。随着计算机技术的快速发展，利用运动目标检测技术对太阳图像进行检测，从而有效地监测太阳活动是目前天文图像处理的研究热点之一。由于受到大气吸收的影响，目前在地面可以观测到太阳的光球层和色球层，而日冕层的观测还仅限于日面边缘。Hα图像可以充分显示色球丰富结构和各类色球活动现象，是最有效的色球观测谱线，国内的观测设备如中国科学院云南天文台全日面Hα色球望远镜[1]、怀柔太阳观测基地的全日面Hα望远镜[2-3]以及光学和红外太阳爆发监测望远镜[4]。国外的观测设备如美国国立太阳天文台的太阳全球振荡监测网(Global Oscillation Network Group，GONG)[5-6]、全球高分辨率Hα网[7](Global High-Resolution H-alpha Network)都使用该谱线观测色球活动。检测太阳Hα图像的特征也就是检测太阳耀斑、日珥、黑子和暗条等太阳活动。

近年来，学者们对Hα图像进行实验并开发了一系列的检测算法和程序。针对太阳暗条的检测方法主要可以分为两类。(1)是基于图像形态处理方法。Hao等[8]采用形态学操作的方法，针对Mauna Loa Solar天文台观测到的Hα全日面太阳图像进行实验，识别暗条的结果较好。Fuller[9]采用了经典的区域生长方法对法国Meudon天文台的Hα图像进行验证，并得到了较好的结果，但是图像局部差异容易导致算法的不稳定性。Gao等[10]采用全局阈值和区域生长方法对美国Big Bear Solar天文台的Hα图像进行验证，此法简单高效，但检测不到边缘暗条，且易漏掉形状较小暗条，黑子也极有可能混杂其中。(2)基于统计机器学习的方法。Zharkova[11-12]利用神经网络方法来识别暗条，并获得了较好的效果。Qu[13]在文献中采用支持向量机的方法来区分暗条和黑子。太阳耀斑的检测也取得了较快的发展，Qu等[14]在2003年利用支持向量机分类器的方法对耀斑进行识别，取得了良好的效果。Borda等[15]利用神经网络算法实现了耀斑的自动检测。Pötzi等[16]针对Kanzelhöhe天文台的Hα图像做临边昏暗改正，并把存在的太阳活动做成高斯直方图，这样就有个灰度从而得到耀斑。从上述可以看出，如何分析海量的观测数据，并开发更有效地算法来检测Hα图像是当今的研究热点。

本文将改进的混合高斯模型方法引入到太阳活动的检测中。预处理图像后，利用三帧差分法先得到的一幅背景图像，再将该图像更新到混合高斯模型的背景模型中，对正在爆发的太阳活动进行检测并识别。通过对观测图像的处理表明：本文算法同时检测出耀斑和暗条爆发等太阳活动，并计算出耀斑爆发时间和级别。通过对中国科学院国家天文台怀柔太阳观测基地和GONG观测的3组数据进行实验，结果证明本算法和现有的算法相比能有效的检测出太阳耀斑、暗条等活动并计算出耀斑级别和时间。改善了混合高斯模型固有的光照变化检测不佳的问题，提高了检测耀斑、暗条等太阳活动的准确度。

1 三帧差分法

帧间差分法是图像处理过程中最常用的一种方法。该算法是将相邻两帧图像之间做差来获得运动区域。通过相邻图像的差值可以快速检测出相邻图像中目标的运动范围。该算法计算速度快，且对于光照变化不敏感。三帧差分法是取连续3帧图像Ii-1(x,y)、Ii(x,y)和Ii+1(x,y)。分别计算两帧间的差值，如式(1)和式(2)所示[17]，选取合适的阈值T得到b(i,i-1)和b(i+1,i)。再对b(i,i-1)和b(i+1,i)做逻辑“与”运算得到运动目标Bi，如式(3)所示。

(1)

(2)

(3)

2 混合高斯模型方法

混合高斯模型(Gaussian Mixture Models，GMMs)广泛应用于数据分析和模式识别等领域，由Stauffer等人[18]在1999年提出。用多个高斯分布(正态分布)来表示图像中t时刻每个像素点的灰度值Xi,t所构成的模型被称作混合高斯模型，其函数表达式如式(4)所示，其中，P(Xt)表示每个像素点K个高斯分布的总概率；η(Xt,μi,t,∑i,t)表示像素点的概率密度；ωi,t表示第i个高斯分布t时刻的权重；μi,t表示第i个高斯分布t时刻的期望；∑i,t表示协方差矩阵；K表示高斯分布的个数，一般取值为3～5，本文取4。随着K值的增大，模型所表示的场景就越复杂，但计算量也随之增加。

(4)

若Xi,t不匹配，当k

ωi,t=(1-α)ωi,t-1+α(Mi,t)

(5)

μi,t=(1-ρ)μi,t-1+ρXi,t

(6)

(7)

ρ=αη(Xt|μi,t,σi,t)

(8)

(2)混合高斯模型的前景检测。首先将每个像素的K个高斯分布根据ωi,t/σi,t比值大小进行排序。再从K个高斯分布中选取B个高斯分布作为背景模型。其中式(9)给出了B的取值。T表示背景阈值。将像素值Xi,t与B个高斯分布进行匹配，若存在匹配，则该像素点为背景点，否则该像素点被检测为前景点。对图像序列中的每个像素点重复执行以上步骤，完成图像的前景分割过程得到Hα图像中的运动物质。

(9)

3 太阳Hα图像特征检测方法过程

3.1 图像预处理

使用改进混合高斯模型检测太阳活动之前，先对Hα图像进行预处理。预处理包括图像的标准化和图像复原两个步骤。本文首先对Hα图像进行标准差归一化。使图像符合标准高斯分布(即均值为0，标准差为1)，保证图像的强度统一。在Hα图像获取时，受到天气条件等各种随机噪声的影响，导致图像存在信噪比低、图像运动模糊等问题。因此要提高图像质量，必须对其进行滤波处理。维纳滤波的原理是假设图像信号可以近似看成平稳随机过程的前提下，使输入图像和复原图像之间的均方误差达到最小的准则函数来实现图像复原。本实验数据的模糊主要是由于天气原因(云，阴天)和望远镜和场景之间的均匀线性运动产生的，因此再利用MATLAB工具箱[21]中的fspecial函数建模并用维纳滤波方法(deconvwnr函数)来对图像进行平滑处理，恢复运动模糊。

3.2 运动目标提取与检测

预处理结束后，采用改进混合高斯模型进行图像的目标(太阳活动)提取。首先确定高斯分布的个数K。本文分别统计了3组数据中多个位置的背景像素的变化，发现3组数据集中大部分的模拟背景像素分别符合4个高斯分布的叠加，因此本文采用4个高斯分布作为数据集的模拟高斯的个数。接着初始化混合高斯背景模型，经过实验，本文采用图像序列第1幅Hα图像做初始背景，以第1幅图像同位置的像素值作为对应的均值，同时方差取比较大的值，第1个高斯分布取较大的权重(本文取0.994)，其他都取较小值来进行初始化。从理论上来说，混合高斯背景建模应该是一种较为完美的背景分割方法。但在发生光线突变时，导致其检测效果下降，为了解决该问题本文将三帧差分法对于光照变化不敏感的优势和混合高斯模型不容易出现“空洞”现象的优势相结合，克服了光照渐变带来的影响。初始化后读入数据，利用三帧差分法先得到运动区域F1，用当前帧与F1做差得到F2，将F2作为背景区域更新到混合高斯模型的主背景分布中建立图像序列中每个时刻图像的动态背景模型[22]。符合高斯分布后进行参数更新。用观测图像减去背景得到前景图像，得到的前景图像就是正在爆发的太阳活动。虽然改善了光照突变带来的影响，但在提取出来的前景图像中依然会存在面积较小的干扰点。这些点会对后续的检测造成一定的影响。由于暗条、耀斑面积通常大于3×108km2[24]，所以小于该值的面积一般为小耀斑，引起显著日地效应的概率不大，所以不是本文的关注对象。(将3×108km2换算到实验数据中图像上像素的大小，经过计算在数据集Ⅰ上面积为376，数据集Ⅱ上为468)。因此本文算法利用MATLAB工具[21]箱中的bwareaopen函数清除掉这些区域，最后完成检测操作。

3.3 太阳活动相关数据计算

(1)判断太阳活动类型。由于暗条的平均亮度要低于太阳日面的平均亮度，而耀斑的平均亮度要高于太阳日面的平均亮度。因此检测出太阳活动后标记在原图像上。计算该区域强度的平均值和整个图像强度的平均值，如果该区域强度平均值小于整个图像的强度平均值，则该区域为暗条活动，反之为耀斑活动；

(2)计算耀斑活动时间。耀斑的开始时刻即是程序识别出耀斑的时刻；峰值时刻即是耀斑面积最大时刻；耀斑面积回到耀斑前状态即为结束时刻。因此利用MATLAB软件工具[21]箱中的regionprops函数分别计算每一帧耀斑(前景图像)的面积来确定耀斑的开始、峰值和结束时刻；

(3)计算耀斑级别。传统上，根据Hα耀斑极大时经投影改正后的面积大小，把耀斑分为4个等级[24]。耀斑等级划分标准[24]如表1所示(R⊙代表太阳半径)。因此将识别的耀斑面积投影改正到日面中心的面积，再换算成平方度对照表1得到耀斑级别。

表1 耀斑等级划分标准

4 实验结果与分析

在Hα的色球观测图像中的太阳活动特征为耀斑、暗条和日珥。耀斑在日面和日面边缘均能发生，暗条是发生在日面上的，在日面边缘的叫日珥，因此挑选出的实验数据集均包含上述特征。本研究选用了来做两个不同观测设备在Hα波段上的观测资料。数据集 I 和 Ⅱ 来自中国科学院国家天文台怀柔太阳观测基地全日面光学和磁场监测系统观测的Hα色球资料。其中怀柔观测的数据所采用的望远镜的主要参数：口径20 cm，有效焦距180 cm，配备半宽0.25 Å的Lyot双折射滤光器，滤光器放在准直光路中，探测器有效像元数为2 712×2 712[23-24]，像元分辨率为0.863 39″。数据集 Ⅲ 是来自Gong的Hα观测资料数据集，其有效像元数为2 048×2 048。数据集 I 包含日面和日面边缘的耀斑和暗条爆发，也包含1个日面上的小耀斑。数据集 Ⅱ 包含日面上的大耀斑。数据集 Ⅲ 包含日面上的中等强度耀斑和暗条爆发。所有数据集的时间分辨率均为60 s。数据集 I 的观测时间是2013/05/12/ 23：01 UT到2013/05/13/ 04：22 UT共255帧图像。数据集 II 是2013/05/14/ 22：42 UT到2013/05/15/ 03：05 UT共135帧图像。数据集 Ⅲ 是从2015/11/04/ 03：10 UT到2015/11/04/ 08：48 UT共309帧图像。

首先对数据集 I 进行实验。预处理结束后设置三帧差分法的阈值T=27。设置GMMs的初始参数：采用4个高斯分布来模拟背景，权值第1位参数取0.994, 其余取等权0.002。方差为25，学习率α为0.005。经过多次实验调整，设置背景分布阈值T为0.7，均值第1位参数取第1帧图像同位置的像素灰度值，其余均为0。经过实验得到结果。图1第1行从左至右分别表示的是数据集 I 中的第146帧(2013/05/13/ 02：26UT)和147帧(2013/05/13/ 02：27UT)的部分原始Hα图像，正在爆发的耀斑活动由圆圈出；图1第2行表示三帧差分法对以上两帧检测出的前景图像；图1第3行表示传统GMMs对以上两帧检测出的前景图像；图1第4行表示本文算法对以上两帧检测出的前景图像。由图1第2行可知，三帧差分得到的结果会有“空洞”现象；由图1第3行可知，虽然检测出了耀斑活动，但会有很多噪声点对计算耀斑数据造成影响，导致识别率降低；由图1第4行可知，本文算法可以较好地检出耀斑区域。

图1 算法比较图

根据检测结果计算得到结论：在日面边缘位置(819″，517″)处编号为NOAA AR 11748的活动区有太阳耀斑爆发。根据耀斑面积随时间的变化计算出耀斑的开始时刻为2013/05/13/ 02：09 UT、峰值时刻为2013/05/13/ 02：32 UT，结束时刻为2013/05/13/ 03：16 UT。经过计算得到数据集Ⅰ的平方度约为9.673，查询表1可知该耀斑级别为2。最后将爆发的耀斑活动位置标记在原始Hα图像上。通过查询Solarmonitor网站，可以得到地球静止环境业务卫星(Geostationary Operational Environmental Satellites，GOES)给出的耀斑数据，如图6第1行。对比两组时间可知，本文的检测结果和GOES卫星给出的时间段大致相同(时间上的偏差由于地面观测时间分辨率设置的不同)。图2中红色部分表示的是耀斑活动(峰值时刻)在Hα图像中的位置。

图2 耀斑活动在原Hα图像中的位置

图3 算法比较图

再对数据集 Ⅱ 进行实验，三帧差分法的阈值T=30。设置GMMs的初始参数：采用4个高斯分布来模拟背景，权值第1位参数取0.994, 其余取等权0.002。方差为20，学习率α为0.006。背景分布阈值T为0.65，均值第1位参数取第1帧图像同位置的像素灰度值，其余均为0。经过实验得到结果。图3第1行左边表示数据集 Ⅱ 中的第73帧(2013/05/15/ 01：51UT)的部分原始Hα图像，正在爆发的耀斑活动由圆圈出；图3第1行右边表示三帧差分法对该帧检测出的前景图像；图3第2行左边表示传统GMMs对该帧检测出的前景图像；图3第2行右边表示本文算法对该帧检测出的前景图像。由图3可知，三帧差分法和传统GMMs会带来识别率降低的问题。本文算法可以很好的避免光照突变带来的误检测现象，并检测出耀斑区域。

根据检测结果计算得到结论：在日面边缘位置(739″，443″)处编号为NOAA AR 11748的活动区有太阳耀斑爆发。根据耀斑面积随时间的变化计算出耀斑的开始时刻为2013/05/15/ 01：36 UT、峰值时刻为2013/05/15/ 01：54 UT，结束时刻为2013/05/15/ 02：16 UT。数据集 Ⅱ 的平方度约为6.521，查询表1可知该耀斑级别为2。查询Solarmonitor网站，得到GOES卫星给出的耀斑数据，如图6第2行。对比两组时间可知，本文的检测结果和GOES卫星给出的时间段大致相同(时间上的偏差由于地面观测时间分辨率设置的不同)。图4中亮点部分表示的是耀斑活动(峰值时刻)在Hα图像中的位置。

为验证该算法的有效性，最后对数据集 Ⅲ 进行相同操作，并得到实验结果。在2015/11/04/ 03：46 UT，日面位置为(167″，34″)处编号为NOAA AR 12443的活动区有剧烈的暗条活动发生。同时伴随着1个较小级别的太阳耀斑爆发，经过计算得到该耀斑的级别为S级。在日面位置(921″，256″)处编号为NOAA AR 12445的活动区有太阳耀斑爆发。耀斑的开始时刻为2015/11/04/ 03：20 UT、峰值时刻为2015/11/04/ 03：24 UT，结束时刻为2015/11/04/ 04：27 UT，耀斑级别为1。图5中亮点部分表示的是耀斑活动(峰值时刻)在Hα图像中的位置，灰条部分表示暗条活动和伴随它发生的小级别耀斑活动在Hα图像中的位置。GOES卫星给出的耀斑数据如图6第3行。本文的检测结果和GOES卫星给出的时间段大致相同。

图4 耀斑活动在原Hα图像中的位置

根据3组数据和GOES卫星的数据对比可以验证本文算法的有效行。表2是3种方法对耀斑检测后计算得到的误检率。根据图1、图3和表2组可知，本文成功检测耀斑和暗条的基础上，改善了混合高斯模型固有的光照变化检测不佳的问题，提高了检测耀斑等太阳活动的准确度。

图5 耀斑和暗条活动在原Hα图像中的位置

图6 GOES卫星给出的耀斑数据(Solarmonitor网站)

表2 3种算法对应的误检率

5 讨论

太阳活动和地球空间环境有着密不可分的关系，对其活动爆发进行识别和跟踪有着重要的意义。本文将改进的混合高斯模型引入到太阳活动检测中，提出一种利用混合高斯背景差分模型建立动态背景来检测Hα图像上正在爆发的太阳活动识别算法。本文分别对3组来自中国科学院国家天文台怀柔太阳观测基地和GONG的数据进行实验。首先对图像进行归一化和维纳滤波处理，再利用三帧差法快速得到背景图像，更新到混合高斯模型主背景分布中建立稳定的Hα图像动态背景，实现运动目标(太阳活动)的提取。本文算法在成功识别太阳耀斑、暗条活动，并计算出了耀斑、暗条活动的日面位置(边缘、日面)。再根据耀斑面积的变化计算出耀斑的级别(C、M、X级)的同时，还提高了检测到准确率。本算法的耀斑检测结果和GOES卫星给出的耀斑级别和时间一致，从而验证了算法的有效性。相比GOES卫星对整个日面的流量监测，本算法不仅可以从Hα序列图像的识别结果中给出耀斑级别，开始、峰值和结束时刻，还可以给出耀斑、暗条等太阳活动特征对应的日面位置，对基于我国现有地基观测的Hα色球数据开展太阳活动实时监测预警具有重要意义。此外，GOES主要是反应太阳日冕活动现象，相比较而言，Hα图像反应的则是更低层次即色球的物理过程，能更好地反应太阳活动的源区特性，因此对于太阳物理有着重要意义。虽然改进混合高斯模型算法的会导致计算量变大，运算时间变长，但是却提高了探测耀斑的准确度。对于耀斑而言提高其检测的准确度比检测速度对监测预警来说有着更重要的作用。