太阳射电II型爆发中频率分裂现象的特征自动提取研究

李国良，袁国武，高冠男

(1.中国林业科学研究院资源昆虫研究所,云南昆明，650233;2.云南大学信息学院,云南昆明，650091;3.中国科学院云南天文台,云南昆明，650011)

0 引言

随着太阳观测仪器的快速发展，大量的观测数据被采集，但是反映太阳爆发活动的数据只占所采集数据的1%[1]。对于观测数据，天文学家通常对太阳射电爆发的部分更感兴趣，而在海量的观测数据中，找出太阳射电爆发就成了一项繁重、枯燥的工作，如果再进一步对爆发中包含的精细结构进行人工提取，就更加重了天文学家的负担。因此，对太阳射电频谱仪采集的数据进行自动射电爆发检测和精细结构的特征自动提取，对于太阳射电天文研究非常有意义。

太阳射电频谱仪对太阳射电进行观测，得到的数据是强度随辐射频率和时间变化的图谱，称为射电动态频谱图。根据动态频谱上的射电信号变化，我们可以找到频率随时间缓慢变化的II型太阳射电暴。1947年，Payne-Scott[2]等人在200MHz、100MHz和60MHz频率上发现了这种具有慢速频率漂移的太阳射电爆发。Wild和McCready[3]于1950年将其定义为II型射电暴。经过多年观测发现II型射电暴通常都展现出基频和谐频辐射结构[4]，在20%的II型爆发中还发现基频或二次谐频会分裂成形态相似的两个部分这一频率分裂现象[5]，该现象的特征参数对于理解II型射电暴的产生机制十分重要。近十几年来，国内外学者仅对少数精细结构开展了特征提取研究，尚未涉及频率分裂现象。所以，目前开展频率分裂现象的特征自动提取研究将会对后续的太阳物理研究工作提供很大帮助。

1 本文采用的太阳射电频谱数据

中国科学院云南天文台拥有两台具有了极高频率分辨率和时间分辨率的全数字化太阳射电频谱仪，能够对更为精细和快速变化的射电爆发及精细结构进行观测。本文实验所采用的太阳射电频谱数据来自云南天文台抚仙湖太阳观测基地的11米口径米波太阳射电频谱仪，这台米波射电频谱仪于2011年7月开始运行，工作频段为70-700MHz，最高频率分辨率为200kHz，优于德国AIP的几百KHz和澳大利亚Culgoora的1MHz分辨率[6]，处于国际领先水平。考虑到时间分辨率和频率分辨率之间的矛盾，该仪器日常工作的时间分辨率被设置为80ms，频率分辨率设置为200kHz。工作时间为北京时间8:00-18:00。

表1 云南天文台米波太阳射电频谱仪的观测参数

2 II型射电爆发区域的确定

2.1 太阳射电频谱图像去噪

本文对原始频谱中的通道做归一化处理，消除强背景干扰和去除通道差异产生的噪声[7]，增强图像中有用的信息。通道归一化的原理是用频谱图像中每个通道的净太阳射电爆发流量来除以该通道的平均值，以此减小各个通道之间的差异来实现噪声消除。由于单通道的净太阳射电爆发流量∆fy(x)为原始通道的流量fy(x)减去该通道的平均值，所以单通道的归一化公式为通过对y取全通道信息实现整幅图像的通道归一化。

2.2 大津法二值化

大津法即最大类间方差法，能够根据太阳射电频谱图像本身的灰度特性,将图像分割为目标前景与背景两个部分。算法公式为分别为前景点数和背景点数所占图像的比例，µa和µb分别为前景和背景的平均灰度。算法能自动选取最佳阈值T计算目标前景与背景的最大类间总方差(T)，实现二值分割。

2.3 形态学重建确定爆发区域

3 频率分裂现象的特征自动提取

经典的精细结构特征提取方法有连片搜索算法[9]、level set提取法[10]、小波变换分离法[11]，这些方法仅适用于区域型精细结构的提取。针对频率分裂现象这类边缘型的精细结构，本文提出了新的方法。

图1 爆发区域确定实验

3.1 高斯卷积平滑图像

高斯滤波是一种线性平滑滤波，常用于图像的平滑处理。本文使用一个标准差为0.5，大小为5×5的卷积核对爆发区域进行二维高斯卷积平滑处理。对频谱图像进行平滑处理可以降低爆发区域噪声，使爆发边界更加清晰，增强精细结构的特征。

3.2 爆发区域的边缘检测

为了消除不良干扰,保留爆发区域的重要结构信息，需要对平滑处理后的爆发进行边缘检测。文献[12]中，P.J.Zhang等人利用Canny边缘检测算法对来自NDA的射电数据进行边缘检测，得到很好的实验效果。但是采用Canny算法对云南天文台的数据进行处理会加强噪声，而且检测出的边缘呈不闭合状态，其他传统边缘检测方法的效果同样不佳，因此本文根据边缘指示函数[13]设计了适合的算法。

3.3 基于距离正则化水平集方法自动提取频率分裂现象

距离正则化水平集是以梯度流形式进行演化的，能够有效处理分裂与合并这一类型的拓扑变化[14]，本文基于Li[14]等人提出的DRLSE算法进行编程实现并应用在基于边缘的主动轮廓模型上，实现距离正则化水平集方法自动提取精细结构。针对频率分裂现象的不规则性，该方法能很好的实现轮廓特征提取，而且演化过程非常稳定。

（1）初始轮廓设置

太阳射电II型爆发发生在频率较低的区域，起始频率不会超过几百兆赫兹。2007年White S针对观测到的频率分裂高频事件报道了频率分裂现象的基频起始频率为380MHz，二次谐频的起始频率为760MHz[6]这一研究成果，这为频率分裂现象的观测和进一步研究提供了一大依据。因此，结合本文实验数据自身的频段信息，将初始轮廓位置设置在380MHz到700MHz的范围内。

（2）频率分裂现象的自动提取

距离正则化水平集方法依赖于图像的梯度和能量函数实现演化功能，能量函数最小化运算驱动轮廓以梯度流形式进行演化，图2（c）中的等高线即为能量值变化的体现，初始轮廓会向着梯度较大的边缘不断进行能量最小化运算以接近边缘。当轮廓越来越接近边缘时，算法会减缓轮廓的演化直至停止。本文根据经验将演化的迭代次数设为固定值1200，能够满足云南天文台数据的精细结构提取。

实 验定义能量函数为E(ϕ) =µRp(ϕ) +Eext(ϕ)，其中µRp(ϕ)为正则项，能够保持演化的稳定性。外部能量函数Eext(ϕ)中包含边缘项和区域向 ：Eext(φ) =λLg(φ) +αA(φ)，通过最小化外部能量函数能够实现演化功能。其中，区域项αA(φ)能够确定演化方向，边缘项λLg(φ)取最小值时说明轮廓已经演化到目标边界处，演化停止。

（3）精细结构提取结果及分析

实验得到如图2（d）中所示的提取结果，说明本文方法确实能有效实现频率分裂现象的特征自动提取。由于本文数据中存在破坏了II型射电暴完整性的坏通道，导致轮廓在垂直方向上的演化最终停在射电暴的截断处，但这并不影响最终的提取效果。假设没有坏通道的干扰，轮廓在迭代1200次后同样会在恰当的位置停止，最终提取出频率分裂现象，而不会肆意演化。

图2 频率分裂现象提取实验

精细结构中的大部分特征参数只与其所在位置有关与爆发强度无关[15]。所以在提取出频率分裂现象后，可以通过记录频率分裂现象的频率最大值点和最小值点所在的通道和时间来计算精细结构的持续时间和带宽，然后利用线性回归方法计算出频率漂移率[16]。

4 总结与展望

本文的不足之处在于满足研究要求的II型爆发数据不多，所以在提取出精细结构之后只是简单描述了可行的特征参数统计分析方法，希望以后在更多数据积累下，实现完善的特征参数统计分析工作。另外，II型射电爆发中还存在Herringbone精细结构，常常出现在频率分裂现象的基频和二次谐频上，在本文的基础上，可以开展相关研究。