基于概率霍夫变换的太阳射电Ⅲ和Ⅱ型暴自动识别及参数提取

沈发新，高冠男，汪 敏

(1. 中国科学院云南天文台，云南 昆明 650216；2. 中国科学院大学，北京 100049；3. 中国科学院天体结构与演化重点实验室，云南 昆明 650216)

作为一个非常活跃的天体，太阳上经常会有剧烈活动或爆发现象(例如耀斑和日冕物质抛射)。耀斑和日冕物质抛射过程中产生的高能粒子和激波等在日地空间传播，有可能到达地球，由此引发一系列日地物理效应(例如太阳高能粒子事件、地磁暴等)，对人类的航天工程、通讯设施以及长距离输油管道等都有非常重要的影响。因此，为了对空间天气进行有效的预报和预警，减少灾害性空间天气对人类生产生活的破坏，深入细致地研究耀斑和日冕物质抛射过程中产生的高能粒子和激波非常必要。

一般来说，天体射电信号十分微弱，但太阳射电暴的信号略强[1]。太阳射电爆发是指太阳活动期间在射电波段出现的剧烈且短促的流量增强现象，根据在射电动态频谱上展现的不同形态，可以分为太阳射电I-V型射电暴。太阳射电Ⅲ型和Ⅱ型暴是目前研究比较多的两种类型，分别是由耀斑和日冕物质抛射期间产生的高能电子束流和激波引起的。

具体来说，耀斑和日冕物质抛射期间的磁重联加速电子形成高能电子，高能电子束流以亚相对论速度(约为0.2～0.6倍光速)沿磁力线快速运动引起周围等离子振荡产生射电Ⅲ型暴，因此，射电Ⅲ型暴在射电动态频谱上表现出快速的频率漂移，是高能电子束流的最佳示踪器。

日冕物质抛射在日冕和行星际空间快速运动时，当日冕物质抛射的速度超过本地的阿尔芬速度时，会产生日冕激波或行星际激波。通常认为太阳射电Ⅱ型暴是由于激波在日冕中向外传播引起等离子体振荡产生朗缪波，再转换为电磁波辐射并以本征等离子体频率和二倍频向外辐射。所以，射电Ⅱ型暴的观测特征一般具有基频和二次谐频结构，由于激波的运动速度远低于高能电子束流，因此，太阳射电Ⅱ型暴在射电动态频谱上通常表现为缓慢的频率漂移。射电Ⅱ型暴可以视作日冕物质抛射激波的最佳示踪器。

目前，通常认为太阳射电Ⅲ型和Ⅱ型暴的主要辐射机制为等离子体辐射，根据等离子体辐射的特点可以通过射电暴的信号频率直接得到射电辐射源所在区域的密度，再根据日冕密度模型可以获得辐射源所处的日冕高度，进而获得辐射源的运动速度，例如日冕物质抛射激波速度[1]。

现有的太阳射电观测设备普遍有较高的频率分辨率和时间分辨率，数据量较大，如果仅采用人工识别的方法，无法满足实际需求。到目前为止, 世界上已经发展了多种太阳射电暴识别算法。这些算法从不同角度处理观测数据，如对流量数据使用阈值法确认有无爆发[2]，对频谱图像数据使用霍夫变换和主动轮廓的方法确定爆发事件的频率-时间关系[3]。

本文在这些识别算法的基础上，使用概率霍夫变换的方法提取太阳射电Ⅲ型和Ⅱ型暴的起止时间、起止频率、频率漂移率等重要参数，提高了识别的准确率和速度，并利用识别出的射电暴的物理参数结合相应的物理模型估计了日冕物质抛射激波的速度。这不但对空间天气的预报和预警有重要的实用价值，而且对研究高能粒子束流和日冕物质抛射在日冕和行星际空间传播物理过程也具有重要的意义。

1 射电Ⅲ型暴的自动识别

射电Ⅲ型暴是一种最常见的太阳射电暴，是太阳剧烈活动后产生的非热高能电子束在日冕和行星际空间运动的最佳示踪器[4]，同时也能最先反映太阳暴能量释放和粒子加速等物理过程，是耀斑、日冕物质抛射等太阳爆发事件的先兆。

由图1可知，在不同的频率范围内，太阳射电Ⅲ型暴的频谱特征不同。例如，低于10 MHz，由于Ⅲ型射电暴辐射源所处的日冕密度随着日冕高度的变化较慢，因此，10 MHz以下的射电Ⅲ型暴的频率漂移速度远低于10 MHz以上的。而高于10 MHz的射电Ⅲ型暴具有类似直线的结构，这使得本文使用概率霍夫变换算法进行射电Ⅲ型暴识别成为可能。

图1 Phoenix 4射电频谱仪、南希10米波射电频谱仪和WIND/WAVES射电探测器在2014年2月16日获得的一例超宽带射电Ⅲ型暴的综合频谱图[5]

本文使用法国南希10米波天线阵(Nancay Decameter Array, NDA)的太阳射电观测数据进行Ⅲ型暴识别，其时间分辨率为1 s，观测频率范围为10～80 MHz，其中均匀分布着400个通道。

为了识别观测数据中包含的太阳射电暴，我们需要对观测数据进行处理。处理过程 (图2) 分为3部分：(1)原始动态频谱图分段及背景去除；(2)动态频谱图前景提取及噪声处理；(3)识别Ⅲ型暴事件并提取爆发信息。

图2 Ⅲ型暴识别流程

1.1 太阳射电频谱数据预处理

在太阳射电动态频谱图上，Ⅲ型暴表现为具有快速频率漂移的信号，且不同Ⅲ型暴事件的相对强度有很大差异，而背景信号强度在几小时内逐渐变化。因此，我们把原始数据切分为若干小片段，这样有助于去除背景效应，同时也能尽量避免信号强的射电暴和弱的射电暴在同一个数据段中，导致弱的射电暴不能识别。经过多次试验测试，我们得出每个片段为10 min的观测数据可以有助于识别并提取射电暴的参数。

在南希10米波天线阵可观测的频率范围内，我们将每个数据段的时间长度设置为600 s，采样间隔500 s，即两个相邻数据段的重叠时间为100 s。这种切分方式可以保证在某个片段中恰好切断的爆发事件在下一个片段中可以识别。 因此，切分之后的数据段都是[600 × 250]的二维数组I(ti,fj)。

由于观测设备不同通道之间存在明显的增益差异，由此产生的通道效应在频谱图上表现为高亮的横条纹，这些横条纹造成爆发结构间断，有时甚至遮盖较弱的爆发事件，影响对太阳射电暴的识别，所以必须先消除观测数据中的通道效应。具体处理过程为在各通道的数据上减去该通道所有数据的均值，以去除背景，只保留来自太阳的射电信号，实现各通道数据归一化，具体公式[3]为

(1)

其中，I(ti,fj)和I′(ti,fj)分别表示处理前和处理后的数据；a和b分别表示该数据段的开始时刻和结束时刻的索引。对数据段I中的每个值进行计算，生成消除通道效应后的I′。

1.2 Yen算法二值化及图像腐蚀

对太阳射电动态频谱图进行预处理后，我们还需要对图像进行二值化，将图像灰度值置为0或255，经过归一化后将灰度值置为0或1，使目标前景与背景区分开来，以消除背景的影响，有利于寻找目标前景中的太阳射电Ⅲ型暴事件。

太阳射电动态频谱图二值化算法有很多种，但效果各有优劣。本文使用多种自动确定阈值的算法对观测数据进行测试，生成二值化图像，原始图像与使用各种阈值法进行二值化的图像如图3。

在图3中，Original代表原始图像，其余图像代表各种算法(如Isodata, Li, Mean, Minimum, Otsu, Triangle和Yen)对原始图像进行二值化的输出结果。其中最为典型的大津法(即Otsu法)[6]，也称为最大类间方差法，是一种使用较为广泛的自适应阈值化方法。该算法通过不断尝试各种阈值，计算图像的类间方差，选取类间方差最大时的阈值作为对该图像进行二值化的标准。此时，原图像以该阈值为界限，将图像中的所有像素分为背景像素和目标像素两大类。

图3 各种自动确定阈值算法的二值化效果

从图3的测试结果可以看出，只有Yen算法[7]能较好地确定阈值，在二值化后明确区分图像前景与背景，并较好地保留前景信息，即爆发的信息，未与背景信息混淆。而其他算法的输出结果未能准确区分图像前景与背景，导致爆发信息被干扰，难以识别。

Yen算法是一种自动多门限阈值法，能够根据太阳射电频谱图像的灰度特性，将图像分割为目标前景和背景两部分。与Otsu法类似，Yen算法能自动选取最佳阈值，实现图像的二值分割，但其选择最佳阈值的依据是最大熵准则，基本思想是通过选择阈值使背景和目标共同提供的信息量最大。

如图4，对太阳射电动态频谱图进行二值化后，我们还需要对图像进行腐蚀(Erosion)操作，去除图像中约每小时一次的定标信号，以及某些在频谱图中占面积较小的瞬时干扰信号。

腐蚀是一种基本的形态学运算，可以寻找图像中的极小区域，将图像中的高亮区域进行缩减，在输出的图像中，高亮区域比原图的更小。

从图4可以看出，腐蚀算法较好地去除了频谱图中的定标信号与瞬时干扰。同时，受腐蚀算法的影响，爆发事件结构的边缘也遭到腐蚀，但整体结构与频率漂移率不变，不影响后续识别。

图4 经过二值化与腐蚀后的图像Fig.4 Image after binarization and erosion

1.3 使用概率霍夫变换算法识别太阳射电暴事件

霍夫变换是一种特征检测方法，广泛应用于图像分析、计算机视觉等领域。霍夫变换(图5)可以在图像中寻找某些图形，如直线、圆和椭圆等。

图5 霍夫变换原理[8]Fig.5 Hough transform principle[8]

常用的霍夫变换检测直线的方法是利用

ρ=xcosα+ysinα

(2)

在图像空间和参数空间之间,将图像中的曲线转换为参数空间的一点,建立对偶变换。(2)式中，ρ为极径；α为极角，取0～180°；x为像素点相对图像原点的行坐标；y为像素点相对图像原点的列坐标。

如图6，输入图像从图像空间转换到参数空间后，在参数空间存在3个明显的亮点，而参数空间的亮点代表输入图像中的直线结构。利用霍夫变换的直线检测功能，我们可以在太阳射电频谱图像中寻找拥有类似直线结构的太阳射电暴。

图6 霍夫变换的效果Fig.6 The effect of Hough transform

本文使用概率霍夫变换算法解析曲线的几何特性，尽量避免对图像中的每个非零像素进行计算量大的投票过程，而是从点集中随机选取一个像素点进行分析。

概率霍夫变换算法具体步骤为(1)每个区间对应一个累加器A(r,θ)，从每个区间的点集中随机选取一个像素点，在参数空间计算θ对应的r值，并在对应的累计器A(r,θ)上加1；(2)从点集中删除该点并更新累加器；(3)若累加器的值大于设定的阈值，则表示已检测到直线，然后删除该直线上所有的点；(4)重复以上步骤，直至点集为空。

根据太阳射电Ⅲ型暴的频谱信号特征，爆发事件在太阳射电动态频谱图中近似一条倾斜的直线。直线的倾斜程度与太阳射电Ⅲ型暴的频率漂移率相关。

我们对太阳射电动态频谱图使用概率霍夫变换寻找Ⅲ型暴。如果在动态频谱数据段中存在符合约束条件的直线，则认为该动态频谱数据段中可能存在Ⅲ型暴事件。但是，在数据段中也可能存在其他类型的射电暴，例如太阳射电Ⅱ型暴在频谱上也可能具有类似直线的图像特征。所以需要根据识别的频率漂移率来限制和过滤其他类型的爆发事件。

动态频谱数据段经过二值化、腐蚀处理后，对每一帧二值化后的图像进行概率霍夫变换，检测图像中的直线结构。同时，已知在南希10米波天线阵观测数据的频率范围内，太阳射电Ⅲ型暴的频率漂移率约为1～10 MHz/s，所以我们需要对检出的直线进行筛选，只留下频率漂移率符合Ⅲ型暴的直线。若经过筛选后，数据段中仍存在一条或多条直线，则确定该数据段中至少包含一个事件。

如图7，图中有直线特征的结构被识别并标红，此处很可能是一个Ⅲ型暴事件，同时还可以获取线段两端的坐标值，即爆发事件的起始时间、起始频率、终止时间和终止频率。图8中事件的起始频率为50 MHz，终止频率为27 MHz，持续时间为15 s，频率漂移速率为-1.5 MHz/s，负值表示Ⅲ型暴从高频向低频漂移。

1.4 南希10米波天线阵数据中Ⅲ型暴的识别结果

以上识别算法可以高效地在大于10 MHz的频谱数据中识别形态类似直线的Ⅲ型暴事件，并获得起止频率、起止时间和频率漂移率等信息，为研究爆发事件提供支持。

如图8，图8(a)中爆发事件的起始频率为67 MHz，终止频率为23 MHz，持续时间为17 s，频率漂移率为-2.6 MHz/s。图8(b)中爆发事件的起始频率为37.3 MHz，终止频率为28 MHz，持续时间为6 s，频率漂移率为1.55 MHz/s。以上两例爆发事件的识别结果证明，本文的识别算法可以识别从各频率开始的爆发事件。图8(c)中爆发事件的起始频率为58 MHz，终止频率为46 MHz，持续时间为6 s，频率漂移率为-2.1 MHz/s。图8(d)中爆发事件的起始频率为48 MHz，终止频率为26 MHz，持续时间为18 s，频率漂移率为-1.22 MHz/s。以上两例爆发事件的识别结果证明，本文的识别算法可以识别各种长度的爆发事件。

图8 识别出的不同参数的爆发事件Fig.8 Burst events with different parameters identified

但由于南希10米波天线阵观测的频率范围有限，且在观测范围内的低频频段中存在较强持续不断的干扰，如果爆发事件的终止频率较低，这些干扰信号可能影响算法对爆发事件终止频率的识别。如图9，算法识别的爆发事件的终止频率与干扰信号所在频率紧密结合，这样在识别爆发事件终止频率时有一定的误差。

图9 识别出的低频爆发事件Fig.9 Low-frequency burst event identified

2 射电Ⅱ型暴的自动识别

从太阳射电观测数据中识别太阳射电Ⅱ型暴的方法与识别Ⅲ型暴的方法类似，需要对数据进行如图2的处理过程。但在使用概率霍夫变换算法识别频谱图中的直线结构及起止点坐标后，需要根据Ⅱ型暴的物理特征进行筛选。例如Ⅱ型暴的基波所在的起始频率通常不高(<300 MHz)，频率漂移速率较为缓慢(<1 MHz/s)。

Ne=N0104.32Rs/R，

(3)

其中，N0=4.2×104cm-3；Rs为太阳半径。

基于以上自动提取算法和日冕密度模型，我们对于2012年11月13日云南天文台米波太阳射电频谱仪的观测数据进行处理，识别其中的Ⅱ型暴，提取起止频率、起止时间、频率漂移率等参数，并计算所对应的日冕物质抛射的特征参数(如图10)。

图10 利用观测数据提取日冕物质抛射特征参数Fig.10 Extracting CME feature parameters from observation data

通过概率霍夫变换算法得到该Ⅱ型射电暴的物理参数，并得出对应的日冕物质抛射在02:02:20 UT-02:05:44 UT的速度为870 km/s，而后，LASCO C2日冕仪在02:24:06 UT发现该日冕物质抛射，并获得该日冕物质抛射的速度为851 km/s(数据来源：https://cdaw.gsfc.nasa.gov/CME_list/UNIVERSAL/2012_11/univ2012_11.html)。通过对比发现，这两种方法给出的日冕物质抛射速度非常接近，因此，我们认为利用概率霍夫变换算法自动识别的参数是合理可信的。

3 总结与展望

本文提出了一种基于概率霍夫变换的太阳射电Ⅲ型和Ⅱ型暴识别算法。首先对太阳频谱图像进行预处理，多次试验测试选择合适的时间长度对原始观测数据进行合理切分；对各通道进行均值化，消除通道不均匀的干扰；对图像进行二值化，提取前景信息；对图像进行腐蚀，消除一些窄带射频干扰信号。然后使用概率霍夫变换算法寻找频谱图像中的直线结构，获取Ⅲ型和Ⅱ型暴事件的起止时间、起止频率和频率漂移率等信息，并根据这些信息进一步筛选频率漂移率符合Ⅲ型和Ⅱ型暴特征的直线结构。

本文还利用该方法从2012年11月13日云南天文台的观测数据中自动提取了Ⅱ型暴的物理参数，并结合日冕密度模型获得日冕物质抛射的特征参数。这些参数对空间天气的预报预警具有实用价值。该方法虽然可以识别典型的Ⅱ型暴，但由于Ⅱ型暴容易与其他类型射电暴同时发生，信号有混叠，未来我们将在更多的Ⅱ型暴数据上对该算法进行检验和改进，以达到更好的、对多种数据来源普遍适用的识别效果。

本文结果表明，该方法能有效地识别太阳射电Ⅲ型暴和Ⅱ型暴，获取射电暴的起止频率、起止时间和频率漂移率等信息，并能快速地对长期、大量的观测数据进行识别和分析。同时，概率霍夫变换算法比标准霍夫变换算法的计算量更小，计算过程更短，更适合在观测过程中实时识别太阳射电Ⅲ型和Ⅱ型暴，提取观测参数得出日冕物质抛射激波速度等关键物理信息，实现空间天气预报预警，维护相关人员、设备的安全。