日晕光度计图像日心自动定位方法*

温 骁，刘 煜，宋腾飞，张雪飞，刘念平

(中国科学院国家天文台/云南天文台，云南 昆明 650011)

现代日晕光度计是用于西部太阳选址的重要仪器，它能观测并记录多波段的大气参量。随着观测的进行，观测到的数据量不断庞大，自动处理观测资料也因此成为必要。在对观测数据的处理中，日面中心坐标的确定对数据处理结果具有决定性影响，这是自动处理数据的第一步，也是从观测资料中获取信息最重要的一步。发展了两种方法用于日面中心的自动定位并在文中逐一介绍，最后对比了两种方法所得的日心坐标的差别及其对测量结果的影响。

1 数据介绍

日晕光度计的数据为4个波段，每帧图像大小为765 pixel×510 pixel。在实际采集时一般是每个数据之间间隔15 s，每一分钟采集一套4波段数据。在天气好的时候，一天可观测时间约为10 h，共2 400个数据，大小约1.8 GB。

文中使用的数据是在云南天文台丽江高美古观测站2011年2月25日观测得到的，选取的时间是从8:30到13:30(UT为0:30至5:30)，数据质量非常好，4个波段共1044帧数据。此段时间内天气非常好，天空万里无云。4个波段的观测曝光时间、滤光片通透中心波长、带宽以及透过率列于表1。

表1 日晕光度计的4个波段的观测曝光时间、通透波长、通透带宽以及ND4滤光片在各波段的透过率Table 1 The exposure times，central wavelengths，widths，and throughputs for the four bands of the ND4 filter on the solar-halo photometer

关于日晕光度计仪器参数的详细资料见文［1］，详细数据见文［2］和文［3］。

2 数据处理目的

观测数据中包括日心强度和日晕强度，计算的均是这些参量在CCD上的读数，在下文用日心强度表示日心强度在CCD上的读数，日晕强度表示日晕强度在CCD上的读数。这些数据作为后期处理则需要考虑ND4滤光片的透过率，但这并不是本文考虑的范围，因此并未修正。将实际图像中的日面像记为日面，其中心在图像中的坐标记为日心坐标。

由于数据中含有仪器衍射光，和ND4滤光片支架的影子，这些区域的数据是不需要的，应当去掉;进而分离出日晕区域，用于计算距离日心一定半径处的日晕强度。由于需要求取日心强度以及距离日心一定半径内的日晕强度平均值、日心强度以及定标后的日晕强度(日晕强度与日心强度的比值)，这些都与日心坐标直接相关，如果日心坐标选取错误，将直接导致数据结果出现偏差。因此发展一套自动并且准确的日心坐标寻找方法非常必要。

3 日面中心自动定位方法1:日面总强度法

要找到日面中心，需讨论选取日面以及求取日面中心的标准。

在可见光下太阳极度近似圆形，在理想的观测条件下，日面应该是标准圆形，其中心就是这个圆形的质心(几何中心)。考虑使用日面轮廓的质心为日面坐标，但是在实际观测中有诸多限制:观测条件限制，如系统略微的杂散光，滤光片或CCD上的灰尘，在阳光路径上飘过的尘土，或是看不见的薄云，薄云会影响日面的亮度梯度，进而影响边缘轮廓的选取，ND4之间的多次反射成像。这些都可能使日面非理想成像，造成日心坐标的较大差异。在选取日面轮廓的时候需要自定义一个边界数值，这样加入了人工选择因素。因此使用日面轮廓质心定日心坐标的方法并不可取。

参看图1(a)，太阳日面强度随半径成梯度变化，离太阳中心越远强度越小，在靠近日面边缘处强度变化极度剧烈。

图1 (a)经过日面的一条线上的强度变化。(b)日面区域的等高线。在观测条件好的时候，观测得到的日面非常接近圆形Fig.1 (a)The brightness across the solar disk.(b)Contours on the solar disk.The solar disk in the image is nearly a perfect circle under good condition

如果选取正确的日心坐标和日面半径，则所有在这个半径内的点的强度和会达到最大值，因为只要这个中心点稍微移动，就会使圆环偏离日面，使得内部强度和迅速降低。这一点可以从日面区域的等高线上看出，如图1(b)，靠近日面边缘等高线变化快，说明强度变化快，相比之下日面中心强度变化要慢，图2为日面大小圆环内强度和椭圆环中心变化，偏离实际日心，强度和迅速降低，也说明了这一点。搜寻方法是选取一定范围内的所有点，每个点计算所有距该点为固定半径的圆环内的点的强度和，找到强度和最大的那个点作为日心。这个半径尽量选取日面边缘强度变化最大时对应的值。具体计算步骤如下。

(1)自动寻找日面中心

观测中日面在图像上经常变化。但是日晕光度计自动跟踪系统和人为调整可以保证日面始终处在减光片遮体的衍射光环内。

在程序中选取半径为R=47 pixel，下文讨论选取不同半径的日心识别情况，以及是如何选取这个半径的。

由于日面并未超出减光片遮体的衍射光环，为了减少不必要的计算时间，日心坐标的搜寻就在这个范围内进行。将得到的日面中心记为点(Cx，Cy)。

(2)程序中R选取不同值时的日心识别情况

在不能确定R的取值时，可以先定几个值分别计算日心坐标，后面将有算法用以确定R的取值，取不同R得到的日心坐标见表2。

表2 取不同R得到的日心坐标Table2 The solar-disk center coordinates for different R values

因此选择点(389，263)作为例图最终的日心。

(3)计算离日心一定距离处的日晕强度

图2 自动寻找日心的算法图示每个白色圆环有一个中心，相应的也有在这个圆环内部所有的点的强度总和。图示曲线为经过日面中心的一条线上的点，以及按上述方法求得的强度分布Fig.2 Illustration for the algorithm of calculating the solar central positionEach white circle has a center and a total enclosed intensity.The curve shows the variation of the total enclosed intensity of a circle when its center is moving along the white line through the solar-disk center

将原图去除衍射光以及支架影子部分后，计算数据中所有距日心为r的点强度的算术平均为距离日心r处的日晕强度，记为Isky(r)。去除衍射光和支架影子后的数据见图3(a)，日晕强度随离日心距离的变化曲线见图3(b)。

图3 (a)去除衍射光以及支架影子部分后的数据。(b)Isky(r)随r变化的曲线，曲线中为0的部分是因为去除了衍射光后留下的空隙Fig.3 (a)The image where the system scattered light and the shadow of the three trestles have been removed.(b)Isky(r)versus r.The gap in this curve is due to the removal of the scattered light

(4)日面边缘强度变化最剧烈的半径的确定

使用Isky(r)分析程序中R的选取，在日面边缘处有最大的强度变化。这里做了类似于函数中的一次求导的处理，定义函数:

该函数随r变化曲线见图4(a)，取得最大值时对应的r为程序中所需的R值。此半径确认后，对于处理未变动仪器所取得的数据时均不变。取日心坐标为 (Cx，Cy)和R=47为半径的圆环画在原图上，显示此日心坐标和R的选取很准确。

图4 (a)Isky(r)－Isky(r+1)随r变化的曲线，在r=47处强度变化最大，取R=47。(b)以点(cx，cy)为中心，以r=47为半径，标在原图上Fig.4 (a)Isky(r)－ Isky(r+1)versus r，with the maximum intensity change at r=47.(b)Image with a white circle centered at(cx，cy)and of radius 47

(5)计算机程序完整步骤:

①找到成像系统的中心，去除滤光片支架以及镜筒的散射光;

②找到支架位置，去除其影子在CCD上成像的那部分数据;

③在减光片遮体衍射光环内部，选择半径R=47的圆环内强度和最大值的点为日心，记为(Cx，Cy);

④将去除滤光片支架和镜筒散射光以及支架影子后的数据，计算距日心为r的日晕强度分布Isky(r);计算所有距日心为4.5实际太阳半径和7.8实际太阳半径内的点平均值为Ihalo;计算以点(Cx，Cy)为中心，加上下左右4个点共5个点取算术平均，为日心强度，记为Io。

4 日面中心自动定位方法2:傅里叶变换相关法

该方法原理简要介绍如下，首先定义一个与真实数据大小相同(765 pixel×510 pixel)的数组。然后以ND4滤光片投影中心(x0，y0)(即上述计算机程序第1步的成像系统中心)为初始日心坐标值，采用半径R=47(由方法1中的第4步计算得到)个像素的圆环为日面。数组的日面内部点赋值为1，数组其它像素赋值为0，生成mask数据g(x，y)。其次，类似于方法1，将原图的镜筒口套圈衍射光以及ND4支架影子部分去除(图3a)，使之成为目标数组f(x，y)。最后，利用傅里叶变换中的相关定理，对函数f(x，y)和g(x，y)进行快速傅里叶相关分析，就得到mask数据与目标数据在相关最大时的Δx与Δy值，于是得到日面中心的实际坐标(x0+Δx，y0+Δy)。其中傅里叶变换的相关定理表达式是:

两个程序可以使用C、MATLAB、IDL等语言编写。程序能够得到观测数据中的有效数据信息，经处理后的数据可用于进一步分析。

5 两种方法识别的日心坐标对比

计算了两种方法在相同参数条件下的日心坐标，由于两种方法的判断标准不一样，需要判定识别日心坐标的准确性和两方法识别坐标的差别，以及日心坐标的差别对所获取的数据的影响，对比的指标有日心坐标、日晕强度、日心强度以及定标后的日晕强度。日晕强度和日心强度均按照方法1的算法计算。

对每组数据均对比了以上指标，对比结果见图5～7。

从图5可以看出两种方法求得的日心坐标并不完全一致。为了准确地描述两种方法的差异，将方法1得到的坐标减方法2得到的坐标，对坐标的差异进行了统计，以占总数据的百分比列于表3。

从表3可以看出两种方法得到的结果相同的较多，相差1个像素的非常少，没有相差超过1个像素的情况。经仔细检查程序，认为这种差别是由于两种方法的差别造成的。

从图6、7可以看出由于坐标的差别造成数据结果的差别:日心强度的差别均在1%以下，定标后的日晕强度(日晕强度与日心强度的比值)差别多数在0.5%以下，整体未超出1%。总的来说，两种方法得到的坐标差别对结果的影响非常小。

图5 4个波段所有数据的日心坐标识别对比上面为x坐标的差别，下面为y坐标的差别。横轴是方法1得到的坐标，纵轴是方法2得到的坐标Fig.5 Comparison between the solar-disk center coordinates in the four wavelength bands measured with the two methods The top row is for the x coordinate，and the bottom row is for the y coordinate

表3 两种方法计算得的各波段数据的日心坐标x坐标和y坐标的相差情况Table 3 The discrepancies of the x and y coordinates from the two methods(in percentages)

图6 各波段数据用两种方法求得的日心强度差别的百分比Fig.6 Differences between the solar-center intensity values from the two methods for all four wavelength bands(in percentages)

图7 各波段数据用两种方法求得的定标后日晕强度的差别百分比Fig.7 Differences between the calibrated solar-halo intensity values from the two Methods(in percentages)

6 结论

本文采用两种方法对日晕光度计观测数据的日心坐标位置进行自动计算。第1种方法是基于日面总强度流量最大值原理，第2种方法是基于实际和理论日面模式识别原理。两种日心识别方法虽然本质相近，但实际采用的判断标准不一样。利用它们对同一时段内的数据进行计算，发现取得的日心坐标结果在4个波段均相差甚微(一个像素之内)。在实际情况下得到100%正确的日心坐标是很困难的，另一方面，这两种方法得到的日心坐标的微小差别也是在计算方法本身的理论误差范围内，即不大于一个像素。因此认为这两种方法均可提供有效的日心坐标信息，使用者可以根据自己的偏好选择其中的任何一种方法进行日心坐标的计算。

图8 由方法1计算的4波段日面中心强度(左图)以及定标后日晕强度(右图)随观测时间的变化曲线Fig.8 Variations of solar-center intensity(left)and calibrated solar-halo intensity with observations(right)for the four wavelength bands as calculated by Method 1

另外，从对比2个参量(日晕强度以及定标后的日晕强度)的差别情况看，这种坐标间的略微差别对取得的这些归算参量影响极小。因此，这两种方法均可用于实际计算，用它们计算得到的正确日心坐标奠定了进一步数据处理的基础(日心强度随时间变化的曲线和定标后的日晕强度随时间变化曲线如图8)，2个参量随观测时间的变化显示出观测地大气对太阳光的散射情况，通过计算可以验证大气对太阳光的散射规律，并求得观测地大气的气溶胶、尘埃和水汽等大气成分。

［1］刘顺庆，段辑，张雪飞，等.现代日晕光度计多波段测光系统初步测试结果分析 ［J］.天文研究与技术——国家天文台台刊，2012，9(2):168－175.Liu Shunqing，Duan Ji，Zhang Xuefei，et al.Analysis of the Preliminary Measurements Taken by the Multi-color Photometric System of the Sky Brightness Monitor［J］.Astronomical Research＆ Technology——Publications of National Astronomical Observatories of China，2012，9(2):168－175.

［2］Lin Haosheng，Penn Matthew J.The Advanced Technology Solar Telescope site survey Sky Brightness Monitor［J］.The Publications of the Astronomical Society of Pacific，2004，116(821):652－666.

［3］刘念平，刘煜，申远灯，等.日晕测量与日晕光度计外缘杂散光抑制试验 ［J］.天文学报，2011，52(2):160－170.Liu Nianping，Liu Yu，Shen Yuandeng，et al.Measurement of Sky Brightness and Suppression of Scattering in Sky Brightness Monitor［J］.Acta Astronomica Sinica，2011，52(2):160－170.