FITS格式以扩展性好、数据量大等优点成为天文学中普遍采用的文件格式[1]。美国宇航局的太阳动力学天文台（Solar Dynamics Observatory,SDO）已能采集4 096像素 × 4 096像素、大小为32 MB的原始FITS图像。图像越清晰所提供的信息越多，但同时增加了图像存储的时空开销，天文图像大数量级存储正成为一个棘手的问题。文[2-3]提出分区域提取特征的方法，并得到广泛的应用。许多图像分割方法的提出依然没有解决现有的问题，因为分割结果保留不相关的背景区域，数据集大小没有显著变化，这些将增加后期特征提取的时空开销。
本文尝试对FITS图片集进行一定可行的删减，包括对图片中目标区域的分离，背景区域的删除，这种处理可以在很大程度上减小数据集，减小存储开销，同时太阳活动区域被完好地保留。
受多种因素的影响，导致采集的FITS图像存在多种噪声，包括外部噪声和内部噪声[4]。本文提出基于方形网格结构的太阳现象目标检测方法，在多种噪声干扰的情况下依然能有效标注完整的太阳现象区域，实验结果表明，该方法对噪声具有良好的抗干扰性能。
1 相关工作近年，随着天文数据的急剧增加，太阳特征自动检测和识别的方法得到了很大的发展。
种子区域生长法（Seeded Region Growing,SRG）是区域分割的常用方法，最早由文[5]提出。该方法能将具有相同特征的联通区域分割出来。即使是早期最简单的种子区域生长法也会面临3个难点：种子点的位置选择，种子数目的选取，以及区域生长策略的使用。文[6]结合边界提取提出的区域生长方法将太阳活动按照磁性显著性的不同定义成不同的封闭轮廓，该方法的初始种子位置选自最大磁性强度的像素点。文[7]将被标注区域的质心作为种子结点，进而使用多种子的区域生长方法检测太阳活动区。另外的一些研究[8-9]也提出了结合形态学方法的区域生长算法。
分水岭算法是边缘分割的常用方法，计算方法由文[10]提出。分水岭算法对于区域边缘具有良好的识别能力，并将不同边缘特征的联通区域分割出来。文[11]将形态学重构方法与分水岭算法结合，对粒状结构明显的太阳现象起到较好的分割效果。文[12]提出一种基于多阈值的多级跟踪算法，达到既定阈值水平的像素逐渐加入到相邻的颗粒中，同时保持不同颗粒间的最小距离。这种方法是以分水岭算法为基础实现的。文[13]提出另一种标记控制的分水岭方法，用来自动检测粒状结构。
阈值分割是最常用的图像分割方法之一。文[14]提出自适应的暗条活动分割方法，暗条的中心点位置和半径的检测采用级联Hough变换圆。文[15]结合区域生长法和全局阈值分割暗条。文[16]利用全局阈值、局部阈值结合的数学形态学方法分割暗条活动。
种子区域生长法和分水岭方法因各自特点被广泛用于太阳活动分割，尤其是活动区[6-7]、暗条[8]，这些太阳活动具有明显的灰度特征。然而，日冕喷射、暗条激活等太阳活动，其灰度特征相对不明显，种子区域生长法和分水岭方法对以上太阳活动的分割会造成过度分割或图像空洞两种极端的结果。相比而言，阈值分割方法更简单高效，并且适用性良好。
本文提出一种自动检测太阳活动的方法，该方法在分割太阳活动之后，能将活动区域保存为新的数据集。首先，该方法构造了一种基于方形网格的存储单元用于存储图像的灰度值，然后利用阈值选取策略从阈值区间选择合适的阈值。最后，利用基于方形网格结构的太阳活动目标检测方法实现不同阈值的分割方法分割不同的太阳活动。
2 太阳活动检测方法本部分详细介绍基于方形网格结构的太阳活动目标检测方法的内容，包括方形网格结构以及多阈值方法。
2.1 图像初始化不可见的背景区域（灰度区间0~99）、可见区域中的较暗区域（灰度区间100~150）与较亮区域（灰度区间151~255）共同构成了可见的目标区域。图像初始化处理如下：
|
原始FITS图像灰度范围是0~65 535，将灰度等级从16位降到8位，不仅能减少存储空间，而且能减少图像处理的时间。
2.2 网格结构未经处理的数据集图像的尺寸很大，如果对图像每个像素点进行处理，会面临巨大的时空开销。本文提出的网格划分能在一定程度上缩减这种开销。
步骤1: 网格划分
[ M,N ]表示原始图像大小为 M行、N列。m × n 则代表每个网格图像的大小，给出 m、n 的求解公式：
|
(1) |
其中， INT函数判断数字是否为整形；
步骤2: 图像存储
(1) 建立直角坐标系
水平方向以方形网格的列数为单位，垂直方向以方形网格的行数为单位，对图像的所有方形网格的存储和计算都建立在这个坐标系的基础上，这样对图片的处理将以网格为单位。
[pixel_row,pixel_col] 表示像素点的坐标； [cell_row,cell_col] 代表网格图像的坐标。因此原始FITS图像能够存储在一个四维数组中，该数组定义如下：
|
(2) |
其中， pixel_row[1,m]； pixel_col[1,n]； 
计算复杂性理论强调在运算效率方面，低维的数据结构要比高维的数据结构表现更优[17]。因此，采用三维数组存储原始FITS图像。
|
(3) |
其中， 
如图 1(a)，原始FITS图像被网格划分。图 1(b)表示四维数组结构存储原始图像。图 1(c)表示三维数组结构存储原始图像。考虑到实验采用的数据集图像有1 024像素 × 1 024像素和768像素 ×768像素两种，可分别分割成32行 × 32列和24行 × 24列数量的方形网格，划分后的每个单元格都包含同等数量的像素个数。
|
| 图 1 方形网格结构 Figure 1 The structure of the square grid |
将网格作为最小的划分单位，相对于传统以单个像素作为最小划分单位，具有3点优势：
（1） 对图像的处理过程加快；
（2） 针对不同区域的网格选取不同的阈值，提高处理精度；
（3） 有效防止图像噪声的干扰。
2.3 最小外接矩形最小外接矩形（Minimum Bounding Rectangle,MBR）是图像分类和识别研究中最常用的方法，被用来提取目标区域的几何特征[18]。这里介绍一种基于方形网格结构的最小外接矩形，求解过程如下：
步骤1: 缝隙码
第2.2节中，建立了网格图像直角坐标系，因此，缝隙码定义如下：
|
其中， Pk∈｛P0,P1,P2 … Pn-1｝ 存储每个凸点的坐标值； P0=(u0,v0) 表示起始点的坐标； dl 表示 Pk 移动的方向， dl∈｛0,1,2,3｝ 。
如图 2，目标区域的边界缝隙码表示为
|
|
| 图 2 目标区域边界缝隙码 Figure 2 The crack chain code of the target region boundary |
步骤2: 求解最小外接矩形
dl∈｛1,3｝ 代表垂直方向； dl∈｛0,2｝ 代表水平方向， P0、P7、P14、P21 是缝隙码中的边界极值。最后可求出最小外接矩形的4个坐标点：
|
将上述基于方形网格结构而提出的改进方法应用于以下算法。
输入：一张来自数据集的原始FITS图像
输出：一张包含完整太阳活动的目标区域图像
|
图像归一化的目的是将原FITS图像的像素灰度范围从原来的(0,65 535)限定为(0,255)，从而减小时间开销和计算难度。参考每张FITS图像的灰度特征选取阈值 T1和 T2。对阈值 T 的选取要满足先行条件，即判断图像的外接矩形是否存在。而 Threshold_number和Threshold_density 则取自经验值。结果如图 3。
|
| 图 3 GBTD处理过程 Figure 3 The processing procedure of GBTD |
(d)所占原图比例35%，对比(d)与(a)，处理结果极大地去除了背景区域的干扰，接下来的实验从若干方面验证基于方形网格结构的太阳活动目标检测方法的可行性。
3 实验与分析 3.1 数据集本文采用的数据来源于太阳动力学天文台①，包含8种太阳现象，每种200张，共1 600张，如表 1。
①http://www.cs.montana.edu/angryk/SDO/data/
| Solar phenomena | Number | Size |
| Active Region | 200 | 768 × 768 |
| Coronal Jet | 200 | 1 024 × 1 024 |
| Emerging Flux | 200 | 768 × 768 |
| Filament | 200 | 1 024 × 1 024 |
| Filament Activation | 200 | 1 024 × 1 024 |
| Filament Eruption | 200 | 768 × 768 |
| Flare | 200 | 768 × 768 |
| Oscillation | 200 | 1 024 × 1 024 |
通过观察图像，发现了太阳图像的一些特征，具有明显太阳现象的区域灰度值较大，而不可见的背景区域灰度值较小。常用的灰度级数有8、16、32等，而人眼识别程度最高的级数是 8[19]，可行的阈值范围在100~255之间，将灰度间划分成100~150、151~200、201~255，然后采用8位共256个灰度等级，绘制灰度范围0~255的灰度直方图。
借助灰度直方图的统计结果，阈值区间100~150所占比例高于60%，151~200和201~255所占比例低于40%。这种差异在灰度图像中表现为明显的纹理特征，进而帮助我们选取阈值 T1和 T2。
3.3 评估切割精确度下面的实验给出一种衡量切割精确度的方法。选择各种太阳活动的数据集的50张图片，一共有400副图像，然后手动标注出这些图像中的太阳活动区。标记区域可以是完整的一个区域，或者是离散的多个区域的组合。以手动标记区域为标准，分别比较种子区域生长法、基于方形网格结构的太阳活动目标检测方法、分水岭方法的精确度。结果见图 4和图 5。
|
| 图 4 3种方法的切割精度比较 Figure 4 Comparison of cutting accuracy among three methods |
|
| 图 5 3种方法的时间开销比较 Figure 5 Comparison of time consumption among three methods |
图 4对比了3种方法的切割精度，每个方法包含400组数据，观察这些离散数据点发现，分水岭方法达到了较高的切割精度，其次是基于方形网格结构的太阳活动目标检测方法，最差的是种子区域生长法。实际上，每种方法的平均切割精度，种子区域生长法是0.53，基于方形网格结构的太阳活动目标检测方法是0.56，分水岭法是0.61。 图 5比较了3种方法的时间开销，明显的区别在于基于方形网格结构的太阳活动目标检测方法的时间成本最小，其次是种子区域生长法，最差的是分水岭方法。另外，每种方法的平均时间消耗，基于方形网格结构的太阳活动目标检测方法是0.41 s，种子区域生长法为2.10 s，分水岭方法是2.53 s。
3.4 评估抗噪性图片噪声作为随机误差，是图像处理中最常见的干扰因素。为了测试基于方形网格结构的太阳活动目标检测方法对多种噪声的抗干扰性，选取密度均为0.05的椒盐噪声和高斯噪声分别加入到图像中，分别与未加噪声的原图进行分割相似度的比较（分割相似度是指对噪声图像与原图分割后的图像的相似度），结果见图 6。
|
| 图 6 （a）高斯噪声图片与原图的分割相似度；（b）椒盐噪声图片与原图的分割相似度 |
观察图 6，加入两种噪声后的图片与原图对比，均达到很高的分割相似度，这一结果在图 6（a）中维持在90%及以上，而在图 6（b）中超过95%。出现相似度较低的情况，其原因是多方面的，既跟高斯噪声的分布相关，也与图片中太阳活动的分布特点相关，可通过缩小网格的大小改进实验结果。
4 总结和展望本文根据数据集图像的灰度分布特点，提出一种基于方形网格结构的太阳现象目标检测方法，其中，多种阈值选取策略能应用于不同的数据集，具有一定的可移植性。基于方形网格结构的太阳活动目标检测方法不仅能识别完整的太阳现象，而且有效降低了数据集存储的开销。在接下来的研究中会考虑采用自动阈值选取策略来改进基于方形网格结构的太阳活动目标检测方法，以期获得更优的分割结果。
| [1] | 崔辰州, 李文, 于策, 等. FITS数据文件的检索和访问[J]. 天文研究与技术，国家天文台台刊 , 2008 , 5 (2) : 116 –123 Gui Chenzhou, Li Wen, Yu Ce, et al. Search and location of FITS data files[J]. Astronomical Research&Technology，Publications of National Astronomical Observatories of China , 2008 , 5 (2) : 116 –123. |
| [2] | Banda J M, Angryk R A. Selection of image parameters as the first step towards creating a CBIR system for the Solar Dynamics Observatory[C]//DICTA'10 Proceedings of the 2010 International Conference on Digital Image Computing:Techniques and Applications. 2010:528-534. |
| [3] | Banda J M, Liu C, Angryk R A. Region-based querying of solar data using descriptor signatures[C]//IEEE International Conference on Data Mining Workshops. 2013:1-7. |
| [4] | Burger H C, Sch lkopf B, Harmeling S. Removing noise from astronomical images using a pixel-specific noise model[C]//2011 IEEE International Conference on Computational Photography (ICCP). 2011:1-8. |
| [5] | Yokoya N, Levine M D. Range image segmentation based on differential geometry:a hybrid approach[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence , 1997 , 11 (6) : 643 –649. |
| [6] | James McAteer R. T, Gallagher P T, Ireland J, et al. Automated boundary-extraction and region-growing techniques applied to solar magnetograms[J]. Solar Physics , 2005 , 228 (1) : 55 –56. |
| [7] | Benkhalil A, Zharkova V V, Zharkov S, et al. Active region detection and verification with the solar feature catalogue[J]. Solar Physics , 2006 , 235 (1) : 87 –106. |
| [8] | Gao Jianlin, Zhou Mengchu, Wang Haimin. A threshold and region growing combined method for filament disappearance area detection in solar images[C]//Conference on Information Sciences and Systems. 2001. https://www.researchgate.net/profile/Mengchu_Zhou/publication/242373392_A_Threshold_and_Region_Growing_Combined_Method_for_Filament_Disappearance_Area_Detection_in_Solar_Images_March_21-23_2001_The_Johns_Hopkins_University_Baltimore_Maryland_USA/links/0deec52cb72d87b80d000000.pdf?inViewer=true&disableCoverPage=true&origin=publication_detail |
| [9] | Parnell C E, DeForest C E, Hagenaar H J, et al. A power-law distribution of solar magnetic fields over more than five decades in flux[J]. The Astrophysical Journal , 2009 , 698 : 75 –82. DOI: 10.1088/0004-637X/698/1/75 |
| [10] | Vincent L, Soille P. Watersheds in digital spaces:an efficient algorithm based on immersion simulations[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence , 1991 , 13 (6) : 583 –598. DOI: 10.1109/34.87344 |
| [11] | Feng Song, Xu Zhi, Deng Linhua, et al. Automatic segmentation of granules of the solar photosphere using morphological reconstruction and watershed transform[C]//6th International Conference on Intelligent Networks and Intelligent Systems (ICINIS). 2013:300-303. |
| [12] | Bovelet B, Wiehr E. A new algorithm for pattern recognition and its application to granulation and limb faculae[J]. Solar Physics , 2011 , 201 (1) : 13 –26. |
| [13] | Xie Z X, Yu D R, Zhang J, et al. Properties of magnetic elements in the quiet Sun using the marker-controlled watershed method[J]. Astronomy and Astrophysics , 2009 , 505 (2) : 801 –810. DOI: 10.1051/0004-6361/200810946 |
| [14] | Yuan Y, Shih F Y, Jing J, et al. Automatic solar filament segmentation and characterization[J]. Solar Physics , 2011 , 272 : 101 –117. DOI: 10.1007/s11207-011-9798-2 |
| [15] | Gao J L, Wang H M, Zhou M C. Development of an automatic filament disappearance detection system[J]. Solar Physics , 2002 , 205 (1) : 93 –103. DOI: 10.1023/A:1013851808367 |
| [16] | Shih F Y, Kowalski A J. Automatic extraction of filaments in Hα solar images[J]. Solar Physics , 2003 , 218 (1) : 99 –122. |
| [17] | Ben-Sasson E, Galesi N. Space complexity of random formula in resolution[C]//IEEE Conference on Computational Complexity. 2001:42-51. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.28.7830 |
| [18] | 卢蓉, 范勇, 陈念年, 等. 一种提取目标图像最小外接矩形的快速算法[J]. 计算机工程 , 2010 , 36 (21) : 178 –180 Lu Rong, Fan Yong, Chen Niannian, et al. Fast algorithm for extracting minimum enclosing rectangle of target image[J]. Computer Engineering , 2010 , 36 (21) : 178 –180. |
| [19] | 谢维信, 秦桉. 人的视觉对灰度级别的分辨能力及视觉内部噪声的研究[J]. 航天医学与医学工程 , 1991 , 4 (1) : 51 –55 Xie Weixin, Qin An. The gray level resolution and intrinsic noise of human vision[J]. Space Medicine&Medical Engineering , 1991 , 4 (1) : 51 –55. |