畸变全日面观测图像的修复

杨云飞，冯松，季凯帆

(昆明理工大学云南省计算机技术应用重点实验室和信息工程与自动化学院，云南 昆明 650500)

畸变全日面观测图像的修复

杨云飞，冯松，季凯帆

(昆明理工大学云南省计算机技术应用重点实验室和信息工程与自动化学院，云南 昆明 650500)

欠佳的平场改正或者不理想的光路会导致全日面观测图像的畸变。针对这一问题，提出了一种解决方案，并给出了详细的算法描述。该算法的核心是扣除一个均匀的径向轮廓曲线，然后在残差图像上进行大尺度中值滤波以获得光滑的改正曲面。通过这一算法，可以极大改善全日面图像的畸变，提高图像质量。最后还对该方法的局限性进行了讨论。

太阳图像；全日面；图像修复；中值滤波

CN53－1189/P ISSN1672－7673

全日面观测是太阳观测中的重要部分，其观测图像一直是太阳物理、日地关系及相关学科研究的基础数据，而且对于监测太阳活动和研究日地关系具有重要的意义，因此国内外众多的地面天文台都在进行不同波段全日面的常规观测。如在Hα波段，全球高分辨Hα网(Global High-Resolution Hα Network，GHN)覆盖9个天文台[1]，全球日震网(Global Oscillation Network Group，GONG)拥有6个站点[2]，我国的国家天文台、云南天文台、紫金山天文台、南京大学等也有各自的观测台站。

通常，全日面观测图像的预处理包括两个基本步骤。(1)平场改正:平场主要用于CCD不均匀性的改正，而且由于全日面观测视场很大，因此平场改正还包括校正望远镜光路的不均匀性，如渐晕、CCD快门效应等等。在全日面观测中，平场测量通常采用KLL方法[3]，即通过连续拍摄多幅位于CCD不同位置的全日面像后计算得到一个平场图像；(2)去除临边昏暗:在可见光波段，太阳呈现明显的临边昏暗现象。这对探测、识别和研究太阳表面的各种活动现象是非常不利的，因此通过扣除处理使全日面像在亮度上变成一个比较平整的状态。常用的做法是通过极坐标转换，然后采用中值平均得到一个临边昏暗轮廓，最后再回到笛卡尔坐标系从观测图像中扣除临边昏暗[4－5]。

图1显示了一个质量较好的全日面Hα观测图像和经过预处理以后得到的平均临边昏暗曲线以及扣除临边昏暗后的图像。

图1 (a)全日面观测图像；(b)平均临边昏暗曲线；(c)扣除临边昏暗后的图像Fig.1 (a)A full-disk solar image.(b)The corresponding average limb-darkening profile. (c)The full-disk solar image after the subtraction of the limb darkening

但是有些时候，扣除平均临边昏暗后的图像并不是均匀的，会显示出一些整体的起伏。图2是从全球高分辨Hα网站[1]获取的3幅扣除临边昏暗后的Hα全日面图像，可以看到有一些大尺度的畸变存在。

图2 3幅扣除临边昏暗后不均匀的全日面像Fig.2 Three full-disk solar images with the limb darkening subtracted.The three images still show non-uniform surface-brightness distributions

随着计算机技术和图像处理技术的发展，有些观测台站也研究了一些全日面图像的修复技术[6－8]，并在此基础上开展下一步的太阳特征自动识别工作[9－11]。本文提出一种算法，旨在改善这些宝贵天文图像数据的质量。从实验结果来看，该算法对畸变全日面图像的修复非常有效，而且具有较强的适应性。

1 全日面图像修复处理方法

对于不同的望远镜和不同的观测时段，畸变不但是变化的，而且没有规律可循。对于宝贵的历史资料，通过补充辅助的测量手段去修正这些畸变基本不可能，因此必须从图像本身的特征考虑修复。这些全日面图像的畸变由于主要来自低频，最直接的修复方法就是采用平滑或者低通滤波的方式获得这些图像的不均匀背景，然后从原始图像上进行扣除。对于全日面图像，无论是色球上的耀斑还是光球上的黑子，它们和日面平均强度有非常大的反差，而且其强度及面积变化也非常剧烈。因此，需要选择一种尽可能和信号强度无关的平滑方式。经过比对实验，发现大尺度的中值滤波是一种较为可行的方案。

还有，由于临边昏暗的存在，必须先将全日面图像变得较为平坦，否则对非日面中心的处理就会极为困难。延用常规扣除临边昏暗的办法展平图像。

全日面图像修复的处理程序包括如下10个步骤:

(1)采用图像边缘检测方法，提取太阳边缘，并使用圆曲线最小二乘法拟合日面中心和太阳半径；

(2)以日面中心为原点，将图像从笛卡尔坐标系转换为极坐标系；

(3)在极坐标系中，将所有半径方向的曲线(也就是各个径向轮廓)按照角度方向取中值，得到一个较为均匀的平均径向轮廓；

(4)对0.7个太阳半径以内的轮廓曲线进行高阶多项式拟合。考虑到太阳轮廓是以日心为对称，而且在日心区域是平滑过渡的，因此多项式的选择不但以日心对称，而且在原点上也是高阶导数连续。标准的临边昏暗曲线较为平稳，但由于要处理的对象是畸变图像，因此选取的多项式高达6阶；

(5)对离开日心0.6个太阳半径以外的径向轮廓进行多点平滑，从而保证整个曲线的光滑性和连续性；

(6)用这个光滑的平均轮廓曲线替代在极坐标平面上的所有径向轮廓曲线，得到一个光滑的曲面，并将其变换回笛卡尔坐标系，得到一个以日心径向对称的全日面像背景轮廓；

(7)按照平场的处理方式，将原始图像除背景轮廓得到一个残差图；

(8)对残差图像进行大尺度中值滤波，从而在残差图像上扣除各种太阳活动特征，得到一个平滑的修正曲面。经过测试，选择太阳直径的2%～10%作为中值滤波宽度较为可行。同时考虑到太阳边缘的情况，在滤波之前，要将图像按照边缘对称方式扩展半个滤波宽度，即便太阳边缘很靠近图像边缘的时候，也可以保证进行滤波处理；

(9)用原始图像除平滑的修正曲面，得到对称的全日面图像；

(10)采用标准扣除临边昏暗的办法获得最终的处理结果。

2 畸变全日面像的处理结果

图3给出了一个畸变的Hα图像经过上述处理的主要过程及结果。从图3的(a)和(b)可以看到，这幅图像是不对称的，有一侧明显发亮，这也反应在极坐标平面图像(c)上；(d)显示了8条间隔45°的太阳径向轮廓曲线，从中可以看到，它们有较大的弥散，这也意味着该全日面像的不对称性；(e)是对所有径向轮廓取中值以后的平均轮廓曲线，可以看到日面中心并不是这个轮廓的最大值，也就是该畸变图像反映的是“临边增亮”现象；(f)是扣除这个轮廓以后的残差图像，其不均匀性明显显现出来；(g)是经过以5%太阳直径为滤波宽度的中值滤波后的光滑残差曲面，也就是畸变改正曲面；(h)是经(g)图修正后的原始图像。可以看到，经过修正的图像呈对称型；(i)是最终的结果，图像表现得非常平坦，并且清晰展现了日面的各种活动情况。这样就可以使用一些简单的阈值法分割暗条、谱斑等特征，而且太阳边缘的日珥也得到了保留。

为了进一步验证本算法的通用性，刻意寻找和选择了一些畸变非常严重的图像，包括平场改正错误和在仪器测试阶段的某些渐晕严重的全日面像。

图3 (a)原始图像；(b)原始图像的三维展示；(c)极坐标图像；(d)45°间隔的太阳径向轮廓；(e)平均径向轮廓；(f)扣除平均径向轮廓后的残差图像；(g)经过大尺度平滑后的修正图像；(h)原始图像除以图像(g)的结果的三维显示；(i)按照标准扣除临边昏暗流程对(h)处理后的结果Fig.3 (a)A raw image；(b)A 3-dimensional representation of the raw image；(c)The raw image under the polar coordinates；(d)Radial profiles within angular intervals of 45 degrees；(e)The average radial profile；(f)The residual of the subtraction of the average radial profile from the raw image；(g)The image after the smoothing of the residual with a largescale median filter；(h)A 3-dimensional representation of the result of the division of the raw image by the smoothed image in (g)；(i)The final result of applying the standard procedure of subtracting a limb-darkening profile from the image in(h)

图4的第1行包括3幅图像，前两幅为Hα像，第3幅为光球Fe I像。第1幅的畸变是由于错误的平场造成的，明显呈现出一种“临边增亮”现象。后两幅图的畸变都是由于光路上的渐晕造成；第2行图像是处理上述3幅图得到的均匀径向轮廓，显然它们与标准的临边昏暗轮廓有较大的不同；第3、4行是扣除这个均匀径向轮廓后的残差图像及其得到的改正图像，大块的不均匀性明显显现出来；最后一行是扣除这些不均匀性和临边昏暗后得到的最终结果，可以看到图像已经变得非常平坦。

图4 不同畸变全日面像的处理过程和结果Fig.4 Intermediate and final results of processing full-disk solar images with different distortion patterns

3 结 论

从上述处理过程的讨论和实际图像的处理结果来看，本文的算法对处理畸变的全日面图像比较有效。不但能处理由于糟糕的平场改正导致的图像降质，还能处理一些光学系统带来的畸变。整个算法的流程并不复杂，易于实现。

然而，该算法的核心是滤波参数的选择，不同的参数会对结果有不同的影响。上述所有图像的处理，均采用了5%太阳直径作为中值滤波器的宽度，也就是可以改正大于这个尺度的畸变，也基本可以保证各种太阳活动不会被平滑掉。但即便这样，在图4的第4行图像上，隐约可以看到在对应的太阳活动部分有些亮斑存在，虽然这些亮斑的起伏已经在1%以下。还有，该滤波尺度对于一些高频变化的畸变显得力不从心。图5是一个快门效应明显的全日面Hα像，在上述滤波尺度上依然可以看到清晰的快门效应，直到将滤波尺度降低到2%太阳直径时，快门效应才消失，但这也意味着这个尺度以上的太阳活动会被当作背景处理。

图5 (a)带有快门效应的全日面图像；(b)5%太阳直径为中值滤波宽度的处理结果；(c)2%太阳直径为中值滤波宽度的处理结果Fig.5 (a)A raw full-disk solar image with the shutter effect；(b)The result of processing the raw image with a median filter of a window width of 5%of the solar angular diameter；(c)The result of processing the raw image with a median filter of a window width of 2%of the solar angular diameter

虽然导致图像畸变的原因有多种，但有些是有规律可循的。比如较为稳定的快门效应和随着望远镜弯沉情况的渐晕，这些都有可能通过对大量数据的处理得到他们系统的变化情况，从而更有利于对这些畸变的修复，这也是未来可以开展的工作。

另外，由于中值滤波是一个非线性滤波，因此处理后的结果并不适合做高精度测量，尤其是高精度的光度测量。而且如果一个活动现象的面积已经达到和图像的畸变相近的尺度时，该图像处理的办法对于低频的信号也会带来一定的影响。但从形态学的角度去判断和识别太阳活动现象，这种改正方法基本能满足需求。

必须强调的是，这些畸变图像是质量较差的图像，而图像处理手段的应用也只能起到尽可能改善图像质量的作用，不可能从根本上恢复原始的图像信息。因此，在观测时把握和提高图像原始质量才是最为根本的关键。

致谢：感谢中国科学院国家天文台、云南天文台对本实验的大力支持，也感谢GONG项目和GHN项目提供实验图像。

[1] Global High-Resolution H-alpha Network[EB/OL].[2013-03-31].http://swrl.njit.edu/ghn＿web/.

[2] Global Oscillation Network Group[EB/OL].[2013-03-31].http://gong.nso.edu/.

[3] Kuhn J R，Lin H，Loranz D.Gain calibrating nonuniform image-array data using only the image data[J].Publications of the Astronomical Society of the Pacific，1991，103:1097－1108.

[4] Zharkova V V，Ipson S S，Zharkov S I，et al.A full-disk image standardisation of the synoptic solar observations at the Meudon Observatory[J].Solar Physics，2003，214(1):89－105.

[5] Denker C，Johannesson A，Marquette W，et al.Synoptic Hα full-disk observations of the sun from Big Bear Solar Observatory-I.instrumentation，image Processing，data products，and first results[J].Solar Physics，1999，184(1):87－102.

[6] Walton S R，Preminger D G.Restoration and photometry of full-disk solar images[J].The Astrophysical Journal，1999，514(2):959－971.

[7] Toner C G，Jefferies S M，Duvall Jr T L，et al.Restoration of long-exposure full-disk solar intensity images[J].The Astrophysical Journal，1997，478(2):817－827.

[8] Zharkova V，Ipson S，Benkhalil A，et al.Feature recognition in solar images[J].Artificial Intelligence Review，2004，23(3):209－266.

[9] Veronig A，Steinegger M，Otruba W，et al.Automatic image segmentation and feature detection in solar full-disk images[C]//A.Wilson Noordwijk，Netherlands.The solar cycle and terrestrial climate，solar and space weather euroconference.Proceedings of the 1st Solar and Space Weather Euroconference.Netherlands:ESA Publications Division，2000，463:455－458.

[10] Toner C G，Jefferies S M.Accurate measurement of the geometry for a full-disk solar image and estimation of the observational point spread function[J].The Astrophysical Journal，1993，415: 852－861.

[11] Zharkov S，Zharkova V，Ipson S，et al.Technique for automated recognition of sunspots on fulldisk solar images[J].EURASIP Journal on Applied Signal Processing，2005:2573－2584.

A New Approach for Correcting Distortions in Full-Disk Solar Images

Yang Yunfei，Feng Song，Ji Kaifan
(Key Laboratory of Applications of Computer Technologies of the Yunnan Province，College of Information Engineering and Automation，University of Science and Technology of Kumming，Kumming 650500，China，Email:flyingwithcloud＠sina.com)

Due to low-quality flatfields used for corrections or the vignetting effects of telescopes，some full-disk solar images show large-scale non-radial distortions.In this paper we propose a processing procedure to fix this problem.The steps are(1)calculating the location of the center of the sun in a raw image，(2) transforming the raw image with the full-disk centered under the Cartesian coordinates to that under the polar coordinates，(3)calculating the median of the pixel values at each radius from the solar center in the raw image，which yields an average radial profile，(4)fitting the radial profile with a high-order polynomial function，(5)removing the radial profile from the raw image resulting in a residual image，(6)creating a model of the distortion pattern for correction by smoothing the residual image in(5)with a large-scale median filter，and(7)producing the final corrected image by a standard procedure.It is apparent that the key factor for the correction is the value of the window width of the median filter.According to our experience a value of the window width between 2%to 5%of the solar angular diameter is suitable for the correction of a full-disk solar image in most cases.We have processed several distorted images and show the results in the paper.This procedure should effectively repair distortions of full-disk solar images and significantly improve image quality. We discuss some drawback of this approach in the paper as well.

Solar image；Full-disk；Image correction；Median filter

P111

A

1672－7673(2014)01－0054－06

2013－03－31；修定日期:2013－04－12

杨云飞，女，讲师.研究方向:计算机应用与天文图像处理.Email:flyingwithcloud＠sina.com