基于快速鲁棒性主成分分析的日冕喷流自动检测方法*

耿成杰，李润鑫，刘 辉，尚振宏

(1. 昆明理工大学信息工程与自动化学院，云南 昆明 650500；2. 昆明理工大学云南省计算机技术应用重点实验室，云南 昆明 650500；3. 中国科学院云南天文台，云南 昆明 650216；4. 昆明理工大学云南省人工智能重点实验室，云南 昆明 650500)

日冕喷流是太阳外层大气(日冕)中普遍存在的太阳瞬变现象[1]，是沿着开放磁力线运动的等离子体流，在不同位置和不同波段观测到的喷流具有不同的形态学和动态学特征[2]。近年来的研究表明，日冕喷流在分析太阳磁场拓扑结构、加热局地日冕以及加速太阳风等方面可能起着十分重要的作用，同时喷流经常伴随着耀斑、暗条和日冕物质抛射(Coronal Mass Ejection, CME)等太阳活动现象[3]。

文[4]通过随机选择20个X射线喷流和极紫外波段喷流分析得出，这20个喷流都起源于迷你暗条爆发，并且暗条的喷发始于喷流基部亮点附近。在喷流喷发过程中存在迷你暗条的上升运动，在暗条到达喷流顶端之前，测量得到的暗条平均长度等于X射线喷流的平均宽度，这与喷流喷发是由微暗条喷发驱动的观点一致。文[5]首次对日冕爆裂喷流和由它引起的同时发生气泡状和喷流状日冕物质抛射的活动现象进行分析，提出了一个解释日冕爆裂喷流的模型，对日冕爆裂喷流的冷热分量和对应动态相关的日冕物质抛射现象进行解释：喷流冷分量是由喷流底部的暗条喷发引起，并进一步导致气泡状日冕物质抛射；喷流状日冕物质抛射与喷流热分量有关，热分量是限制暗条的闭合场与周围开放场磁力线重新联接而产生的向外移动的热等离子体。文[6]使用太阳动力学天文台卫星的大气成像仪、日震学与磁场成像仪(Helioseismic and Magnetic Imager, HMI)的图像数据和1 m新真空太阳望远镜(New Vacuum Solar Telescope, NVST)的Hα线心图像数据对喷流的演化和喷流源区两个相邻暗条细丝间的变化进行分析发现，该双向环状喷流的爆发是由于喷流源区相邻暗条细丝间的磁重联导致，和目前普遍认为的双向喷流的触发机制有很大不同。

随着太阳观测设备的发展，高时空分辨率的日冕图像数据越来越多，通过图像处理技术从海量数据中自动检测日冕喷流，有助于研究人员对喷流的形成过程、触发机制以及喷流对日冕的加热和太阳风加速背后的物理规律的研究，同时也能够为分析喷流与暗条、日冕物质抛射等伴随活动现象之间的内在联系提供便利。

文[7]将索贝尔(Sobel)边缘检测和霍夫(Hough)变换相结合对太阳动力学天文台的大气成像仪观测数据中的日冕喷流进行检测，并将检测的日冕喷流相应的活动持续时间、长宽度和视线速度等参数存入日球物理学事件知识库(Heliophysics Event Knowledgebase, HEK)。霍夫变换将检测对象从目标空间转换到参数空间，避免了在目标空间检测时的分类运算，但是该方法比较依赖原始图像空间中检测对象的曲线表达式。而日冕喷流在形态上复杂多变，边缘曲线的表达式难以确定，造成检测效果不佳。文[8-9]使用帧间差分法对日出(Hinode)卫星X射线望远镜(X-ray Telescope, XRT)观测数据进行检测，并对检测的喷流长度、速度等参数进行计算和统计分析。帧间差分法原理简单，易于实现，但是对喷流运动速度缓慢和强度较弱时的检测结果不理想。文[2]首先对1 m新真空太阳望远镜观测的多个Hα日浪事件进行研究，得到该类事件的观测特性和主观判断结果，然后对观测数据中比较暗的结构(黑子、暗条和针状体等)进行标记，再利用阈值分割方法分离观测数据中的暗结构，排除已经标记的暗结构，将剩余的暗结构和主观判断结果进行比较，最终得到Hα日浪。该方法需要对黑子、暗条等结构进行标记，比较费时。

事实上，日冕喷流的检测可以看作在日冕序列图像中检测突变，或者是从复杂随机变化的背景中检测有规律变化的喷流活动。因此，如何实现复杂变化背景与稀疏前景的分离是检测日冕喷流的关键。快速鲁棒性主成分分析方法在运动目标检测中有广泛的应用。本文用该方法对太阳动力学天文台的大气成像仪观测数据中不同时间段、不同波段、不同观测位置的日冕喷流进行检测。检测结果表明，该方法与传统的帧间差分法相比能够检出更加完整的喷流，检出的喷流区域孔洞较少，比较连续。

1 太阳动力学天文台/大气成像仪实验数据

太阳动力学天文台是美国国家航空航天局(National Aeronautics and Space Administration, NASA)于2010年2月发射的一颗空间天气卫星，搭载的大气成像仪具有较高的时空分辨率，能够为进一步研究太阳爆发活动提供高质量的观测图像数据。本文使用两组不同时间段、不同波段和不同观测位置的日冕喷流实验数据集。数据集Ⅰ：波段为17.1 nm和19.3 nm，视场范围449 × 179像元，像元分辨率为0.6″/pixel，视场中心坐标(205.2″，-253.8″)。17.1 nm波段的观测时间为2014年10月4日10:10:00～10:41:00(UT)，共156帧图像；19.3 nm波段的观测时间为2014年10月4日10:10:08～10:41:08(UT)，共156帧图像，时间分辨率都为12 s。数据集Ⅱ：波段为17.1 nm和19.3 nm，视场范围186 × 313像元，像元分辨率为0.6″/pixel，视场中心坐标(-219″，-366.9″)。17.1 nm波段的观测时间为2014年3月3日19:10:00～20:05:00(UT)，共276帧图像；19.3 nm波段的观测时间为2014年3月3日19:09:56～20:04:56(UT)，共276帧图像，时间分辨率都为12 s。

这两组数据中喷流变化尺度较大，像素强度增加比较明显，同时还存在一些亮度增强的非日冕喷流变化，导致视场范围内的图像强度发生复杂变化。图1(a)显示了数据集Ⅰ在10:25:00(UT)17.1 nm波段对数对比度增强后的全日面图像；图1(b)显示了(a)中矩形区域视场范围内17.1 nm波段在该时刻对应的对数对比度增强后的日冕图像；图1(c)显示了19.3 nm波段在10:24:56(UT)对应的对数对比度增强后的日冕图像。(b)和(c)中的日冕喷流在像素强度上分布不均匀，且(b)中的图像内容更复杂，日冕喷流与图像背景像素强度差异大，易于区别。19.3 nm波段在10:27:20(UT)的日冕喷流强度较弱，与图像背景强度差异小，此时17.1 nm波段中的日冕喷流像素强度与图像背景依然有较大的区别。图2显示了数据集Ⅱ在2014年3月3日19:10:00～20:05:00(UT)时间段17.1 nm/19.3 nm波段的对数对比度增强后的日冕图像，图2(b)中的图像内容比图2(c)中的复杂，但(c)中的日冕喷流在直观上易于辨别。

图1 (a)数据集Ⅰ在10:25:00(UT)时刻17.1 nm波段对数对比度增强后的全日面日冕图像；(b)为(a)中红色矩形框标记的日冕图像；(c)为(a)中红色矩形框在10:24:56(UT)时刻19.3 nm波段标记的日冕图像

图2 (a)数据集Ⅱ在19:31:49(UT)时刻17.1 nm波段对数对比度增强后的全日面日冕图像；(b)为(a)中红色矩形框标记的日冕图像；(c)为(a)中红色矩形框在19:31:56(UT)时刻19.3 nm波段标记的日冕图像

2 快速鲁棒性主成分分析方法

鲁棒性主成分分析[10]方法的数学模型为

(1)

(1)式是一个非凸的优化问题，目前还没有针对该问题有效的解决方案，可以通过松弛技术，将(1)式中矩阵S的l0范数使用l1范数代替，使用核范数‖·‖*代替rank(·)近似地逼近矩阵L的秩，由此可以将(1)式松弛为一个凸优化问题进行求解，数学表达式为

(2)

其中，D为观测矩阵；‖L‖*为矩阵L的核范数；‖S‖1为矩阵S的l1范数；λ(λ> 0)为权重参数，用于调节低秩矩阵和稀疏矩阵的占比。假设观测矩阵中含有观测到的m个数据，每一个数据列向量化为矩阵D∈Rm×n的一列，矩阵D可以分解为低秩矩阵L和稀疏矩阵S，即D=L+S。其中L和S未知，L满足低秩特性，S满足稀疏特性，数值可以为任意大小，这也是主成分分析鲁棒性的体现。

由于鲁棒性主成分分析是通过代价高昂的优化来计算的，在数据量较大时，实时性较低，因此鲁棒性主成分分析快速算法的开发显得尤为重要。文[11]对(2)式的惩罚约束进行改进，得到在视频背景建模/前景检测下比非精确增广拉格朗日乘子法(Inexcat Augmented Lagrangian Multiplier, IALM)快一个数量级的快速鲁棒性主成分分析方法，数学模型为

(3)

由于性能的原因，文中主要考虑rank(L)=t时的情况，

(4)

对(4)式使用交替最小化方式进行求解，具体求解过程为

(5)

(6)

(5)式通过对D-Sk进行部分奇异值分解来计算矩阵L，t为奇异值的个数，当t很小时，算法有非常高的计算效率。(6)式通过对D-Lk+1进行逐像素的软阈值操作[shrink(D-Lk+1,λ)]计算矩阵S。软阈值函数定义为

shrink(x,ε)=sign(x)max{0,|x|-ε}.

(7)

快速鲁棒性主成分分析方法是根据低秩和稀疏分解的思想进行运动物体的检测，基本思路是假设由序列图像组成的观测矩阵D可以分解为低秩背景矩阵L和稀疏前景矩阵S。由于太阳日冕时刻发生复杂随机的变化，所以日冕图像中的内容复杂多变。从数据集Ⅰ和Ⅱ中的序列图像来看，日冕喷流在时间序列上变化尺度大，而序列图像中日冕变化尺度比较小的部分在时间序列上差异不大，且在序列图像中有比较高的比例，因此变化尺度比较小的这部分日冕可以近似认为是低秩、相关的。日冕喷流序列图像数据的这一特点与快速鲁棒性主成分分析的低秩、稀疏分解思想相似，因此，本文尝试将快速鲁棒性主成分分析方法用于日冕喷流检测，并在实验对比分析中验证方法的可行性和优势。

3 日冕喷流检测与结果分析

在使用快速鲁棒性主成分分析方法检测日冕喷流运动之前，我们先对Level 1.0的原始图像数据进行预处理。预处理包括FITS(Flexible Image Transport System)文件标准化、太阳自转校正、曝光时间校正和图像对比度增强。对于FITS文件的标准化和太阳自转校正，太阳软件包中提供了相应的解决方法(aia_prep.pro和drot_map.pro)。FITS文件标准化主要实现文件日面中心对齐，像元分辨率统一为0.6″并消除侧倾角。太阳自转校正主要消除由太阳自转在投影平面造成的位移，减小对日冕喷流演化分析时的影响。FITS头文件包含的图像数据曝光时间大小不一定相同，导致图像序列的亮度差异大，把图像数据除以曝光时间得到单位时间内曝光量的图像数据，然后再使用标准差归一化方法，使图像序列的强度统一。原始FITS文件的图像数据对比度较低，这里使用对数函数对图像进行对比度增强，方便进行图像内容观察和喷流检测。

对图像进行预处理之后，本文使用快速鲁棒性主成分分析方法对序列图像进行喷流检测。我们使用快速鲁棒性主成分分析方法分别对数据集Ⅰ和数据集Ⅱ进行实验，对检测结果进行块匹配和3D滤波(Block-Matching and 3D filtering, BM3D)[12]去噪后，使用文[13]提出的最大类间方差法分割图像。并将二值化后连通域面积小于20的区域剔除，得到最终检测结果如图3、图4。

本文以两组不同时间段、不同波段和不同观测位置的日冕图像序列作为检测对象，用快速鲁棒性主成分分析的检测结果与帧间差分法的结果进行对比分析，说明快速鲁棒性主成分分析方法的可行性和优势。从图3、图4中第2行(使用MATLAB的Image labeler应用程序中的像素标记工具，将有喷流的地方直接标记，标记的地方值为1，其余地方为0，将标记结果作为标签使用)和第3行(快速鲁棒性主成分分析检测结果)图像序列对比，快速鲁棒性主成分分析能够检出喷流的绝大部分，且与帧间差分法(第4行，当前时间帧与前12 s相减)检出的喷流区域相比，空洞较少，且检出的区域比较连续，说明鲁棒性主成分分析方法能够用来检测日冕喷流。

对于数据集Ⅰ，从图3(a)中矩形框标记的日冕喷流的检测结果来看，在10:27:12时刻，快速鲁棒性主成分分析方法的检测结果明显优于差分法的，且检测结果中离散点较少，检出的喷流区域更加完整；图3(b)中矩形框标记的日冕喷流的检测结果显示，快速鲁棒性主成分分析方法能够检出喷流的微小区域，差分法则没能检出。

图3 数据集Ⅰ。(a)17.1 nm波段检测对比结果；(b)19.3 nm波段检测对比结果Fig.3 Detection and comparison results of data Ⅰ. (a) 17.1 nm images; (b) 19.3 nm images

对于数据集Ⅱ，从图4中矩形框标记的日冕喷流的检测结果来看，与图3矩形框标记的日冕喷流的检测结果基本一致。图4(a)中矩形框标记的日冕喷流的检测结果显示，快速鲁棒性主成分分析方法能够检出更多的喷流，检测结果中不存在大量的非喷流区域。图4(b)中矩形框标记的日冕喷流的检测结果显示，快速鲁棒性主成分分析方法与差分法相比能够检出日冕喷流，且日冕喷流检测结果从直观上容易辨认，不存在误检的情况。

从图3、图4日冕喷流的检测结果来看，快速鲁棒性主成分分析方法能够在复杂多变的日冕背景中检出运动缓慢且强度较弱的日冕喷流。从原理上讲，帧间差分法对图像内容单一、背景稳定的图像数据时具有较好的适用性，但是对包含复杂随机变化的日冕图像却存在漏检(图4中矩形框4)的情况。

但是快速鲁棒性主成分分析方法和帧间差分法的检测结果有时相差不大(如图3中(a)和(b)，图4中(a)和(b))，且差分法在检出喷流区域的同时，也检出许多不是喷流的区域。为了对快速鲁棒性主成分分析方法和帧间差分法的检测结果进行客观分析，同时也为快速鲁棒性主成分分析方法中的t和λ的设置提供客观依据，本文采用F1指标[14]对两者的检测效果进行综合评估，公式为

图4 数据集Ⅱ。(a)17.1 nm波段检测对比结果；(b)19.3 nm波段检测对比结果Fig.4 Detection and comparison results of data Ⅱ. (a) 17.1 nm mages; (b) 19.3 nm images

(8)

(9)

其中，P为准确率(Precision)；R为召回率(Re-call)；F1为综合性能指标；TP(True Positives)，FP(False Positives)和FN(False Negatives)分别表示正确检测为日冕喷流的前景像素点、错误检测为日冕喷流的背景像素点和错误检测为背景的前景像素点。综合性能指标F1越高，表明方法的检测效果越好，如表1。对于数据集Ⅰ和数据集Ⅱ，快速鲁棒性主成分分析方法在综合性能指标上均有提高，且日冕喷流检测的准确率明显提高。经过多次实验，根据F1指标高低设置快速鲁棒性主成分分析模型中的t和λ值，数据集Ⅰ中17.1 nm波段对应的t为30，λ为0.01，数据集Ⅱ中17.1 nm波段对应的t为7，λ为0.01。19.3 nm波段采用和17.1 nm波段一样的值。

表1 快速鲁棒性主成分分析与帧间差分检测结果综合比较

此外，在对日冕喷流进行统计研究时，日冕喷流的速度、长度和宽度等参数经常作为统计对象，这里使用快速鲁棒性主成分分析方法的日冕喷流检测结果二值图像的最小面积外接矩形估算喷流的长宽比，继而估计喷流的大致形态。最后根据检测结果的二值图像凸壳质心的位移叠加，估算喷流的总体速度。

根据估算结果得到结论：日冕喷流的长宽比都大于1，形态上具有较明显的朝向性。在数据集Ⅰ 17.1 nm波段测得的速度为174 km·s-1，长宽比为3.19，在10:25:00时达到最大长度约8.5 × 104km，活动时长为10:14:00～10:28:00；在19.3 nm波段测得的速度为194 km·s-1，长宽比为3.25，在19:33:44时达到最大长度约9.7 × 104km，活动时长为10:13:31～10:26:03。在数据集Ⅱ的17.1 nm波段测得的速度为174 km·s-1，长宽比为4.74，在19:33:12时达到最大长度约9.3 × 104km，活动时长为19:25:24～19:34:48；在19.3 nm波段测得的速度为157 km·s-1，长宽比为3.15，在19:32:20时达到最大长度约8.03 × 104km，活动时长为19:22:44～19:34:44。

4 结 论

本文将快速鲁棒性主成分分析方法应用于日冕喷流检测，通过低秩和稀疏分解的思想对日冕的动态背景进行建模，从而检出作为前景变化的日冕喷流。实验结果表明，快速鲁棒性主成分分析方法在喷流变化尺度较大、像素强度增加明显，同时还存在亮度增强的小尺度非喷流变化的日冕序列图像中检测日冕喷流时，检测结果与帧间差分法相比，检出的日冕喷流区域孔洞较少，区域连续，且检测的准确度有所提高。差分法原理简单，易于实现，但是在针对日冕喷流检测这种除了日冕喷流的变化之外，还存在其他小的、亮度增加的非喷流变化的动态场景时效果不理想，而快速鲁棒性主成分分析方法能够检出帧间差分法不能检测到的强度较弱的日冕喷流。

目前，快速鲁棒性主成分分析方法存在的问题是模型中参数t和λ的设置，这些参数还不能实现自动设置，需要通过多次试验人工选择，而如何通过客观指标分析对模型参数进行自动设置和对可行性进行客观评价是改进的重点。日冕喷流的观测特征参数(如速度、长度和宽度等)还需要更加精准的测量。此外，快速鲁棒性主成分分析方法在数据量比较大的情况下运算量大，耗费时间长，如何对全日面上的日冕喷流特别是弱的日冕喷流进行检测是下一步的重点内容。