彭青玉 (暨南大学计算机科学系, 广州 510632; 中国科学院国家天文台, 北京 100012;中国科学院光学天文联合开放实验室昆明基地, 昆明 650011) 摘要 在克服前人图像处理技术缺陷的基础上, 我们成功试验出一种基于土星光环边缘检测的处理技术曀¥精确测定土星位置的方法, 并进一步实施了对土星本体和光环产生的晕的处理. 因此, 可以获得更多的机会来精确测量离土星本体和光环很近的土卫1(Mimas)和土卫2(Enceladus)的位置. 利用中国科学院云南天文台1 m望远镜3个晚上观测的127帧CCD图像实际测量表明: 土星光环几何中心位置的测量与其主要卫星(土卫3~6)位置的测量具有同样好的精度, 单次测定的标准误差优于0.05. 新的图像处理技术将对CCD子午环观测土星位置以及长焦距望远镜观测土星暗卫星位置均有重要的推动作用. 关键词 图像处理 土星光环和卫星 位置测量 土星及其主要卫星的天体测量观测有着悠久的历史, 这种观测的直接目的是改进行星和卫星的轨道理论, 也为当今太阳系空间探测服务. 土星的位置测量曾主要依赖于子午环的直接观测, 但自Sinclair[1,2]独创性地用照相方法通过土星卫星的观测间接测定土星位置以来, 国际上很少发表直接测量土星位置的结果. 主要的原因是土星相位效应的影响, 观测精度不高. Pascu和Schmidt[3]采用了一种以测量土星光环几何中心代替测量土星质心位置的方法,并揭示这一方法与测量土星主要卫星再按卫星轨道理论导出土星本身位置的方法有着同样的 精度(0.20~0.27). 但测量精度仍然不高. 以后, 人们普遍采用以卫星位置测量来导出土星质心位置的方法, 最新的文献有[4~7]. 他们都采用CCD子午环进行观测, 每颗卫星的位置测定 精度(通常为土卫3~6)为0.10~0.20. 以土星光环几何中心位置的测量代替土星质心位置的测量有其优越性. 首先, 由于光环本身很薄, 不存在土星本体那种不断变化的"相位效应"; 其 次, 光环边缘是一个天然的极好的椭圆, 这一椭圆中心与土星本体质心的等价性理论上不存在争议的问题. 但为什么以测量土星光环几何中心来代替土星质心位置的方法没有得到很好 的发展呢 这可能与下述两方面的原因有关: () 高精度测量光环的边缘点是一个棘手的问题, 因为地面观测的光环图像边缘是模糊而非清晰的; () 由于土卫空间探测的需要(例如CASSINI)人们更多地把注意力集中于土星卫星的天体测量. 例如, 一些长焦距望远镜照相或CCD观测投入了使用, 人们开发出了更多的测量土星卫星的技术, 测量精度也在不断提高. 目前长焦距望远镜CCD观测土星主要卫星(土卫3~6)的精度(单次测定的标准误差)约为0.08[8,9], 更好的可达0.04~0.06[10,11]. 应当看到, 无论是改进土星轨道理论或者是改进土星暗卫星(例如土卫1, 土卫2)轨道理论都需要精确测定土星本身的位置. 一方面, 当用CCD子午环等短焦距望远镜观测土星时,由于CCD视场较大(如美国海军天文台FASTT的视场为0.7平方度)能够找到足够的定标星(如ACT星). 如果能测定土星光环几何中心的位置, 则土星的位置将不受卫星理论位置误差的影响; 另一方面, 用长焦距望远镜测定暗卫星位置时, 需要知道土星中心的位置以便能精确扣除行星晕对暗卫星位置的影响. 人们曾研究开发过这种晕的处理技术, 但不都是令人满意的.例如Pascu等人[12]从图像分析方法入手将近主星周围晕的光度分布模拟为二阶多项式以便对卫星所在位置的晕进行扣除. Stone和Harris[4]的处理方法与Pascu等人[12]的方法有些相似性. 他们假定行星晕的分布满足圆对称性(见文献[4]中方程(2)和(3)),而且模拟的多项式高达6阶项. 对离主行星较远的卫星, 这种假定也许能较好地满足, 但对于具有明显椭圆光环结构的土星星像, 这种假定无疑是过分简单的, 特别是考虑近主星范围内的晕时尤其如此. 此外, 高阶拟合项的有效性值得怀疑,因为这容易让人联想到数值计算中的龙格现象. 文献[11]的处理方法虽然克服了上述缺陷, 但采用轴对称方法(见文献[11]中图 2)致使对称轴上的星像的处理有时变得不方便, 而且这种单个像素精度的中心位置, 对近主星附近的卫星所在位置的晕的去除难免担心其引起附加的非对称性. 再者, 该方法去晕后的背景灰度水平(零附近)很低, 人眼难以识别卫星附近的背景, 妨碍了人机交互处理. 本文从土星光环边缘检测的角度入手, 逐步扣除影响精确测定光环几何中心位置的各种因素, 并与最新最佳土星卫星理论TASS1.7[13,14]进行比较以验证测量方法的精确性. 因此, 本文的研究并非提供土星相对于恒星的天体测量位置, 而是研究土星光环几何中心位置精确测定的关键技术. 此外, 我们还实施了精确的去晕技术以及设置去晕后背景光度水准的办法, 以方便卫星位置测量的人机交互处理. 1 试验用观测资料 2002年1月1~4日, 我们在中国科学院云南天文台1 m望远镜上获得了3个晚上(1月2日晚由于天气关系没有观测资料)共计127帧土星及其卫星的CCD图像. 观测过程中均使用I型滤光片. 每晚上取得的观测资料数及CCD露光时间等列于表1. 所用望远镜及CCD器件 的参数可参考文献[11]中的表1. 1月3日晚有两个观测资料子集的原因是观测过程出现了机械故障, 我们调整了安装CCD的底座, 从而定标系统发生了变化. 表1 观测资料说明 资料集观测日期图像帧数露光时间/s 0201011月1日271.5~2.5 020103a1月3日211.5~2.0 020103b1月3日251.8~2.0 0201041月4日542.0 比较表1与文献[11]中表2可以看出, 本次观测采用的露光时间相对较短, 目的在于想获得未过饱和的土星及光环的星像以便于光环图像的边缘检测. 2 土星光环及其晕的图像处理 2.1 土星光环的处理 图1 一帧典型的土星及其卫星的CCD图像 图1是一帧典型的土星及其卫星的CCD图像,从图中可以见到土卫2(Enceladus)土卫3(Tethys)和土卫4(Dione), 而离土星最近的土卫1(Mimas) 淹没在土星及其光环的晕之中. 为了得到精确的土星光环几何中心的位置, 我们采用经典的置阈技术(即设置一个合适的灰度阈值将原始图像转化为一个二值图像), 再实施边缘检测获得光环的边界, 最后采用椭圆拟合方法获得光环几何中心的位置. 具体的执行步骤如下: 第1步 选取一个包围土星及其光环的矩形区域(图2(a)), 搜索矩形区域内的灰度值获得最大最小灰度值M和m. 设阈值V有下面的形式: ()VmMm=+-k,(1) 其中k由试验确定. 用文中脚注1)的方法表明k = 0.5时能获得较好的效果. 第2步 二值化上述矩形区域内的子图像, 并用一个四邻点模板(见文中脚注1)的附录A) 检测获得子图像的边缘(见图2(b)). 我们看到除了土星光环的外面边缘被检测到以外, 也同时检测到了土星本体和光环之间的边缘部分. 对于要求拟合的椭圆而言, 后一部分是多余的. 第3步 改进边缘检测后的二值化图像, 得到光环边缘. 具体地, 我们在水平和垂直两个方向分别进行每一行和列的搜索来剔除内部点. 处理后的图像见图2(c). 第4步 用椭圆方程迭代拟合上述边缘点, 并对偏离拟合曲线上较远的点进一步剔除.具体地, 设椭圆的两焦点坐标为(X1, Y1)和(X2, Y2), 半长轴为a, 则任一边缘点(X, Y)到两焦点的距离为 2222 1122()()()(),DXXYYXXYY=-+-+-+-(2) 我们知道, 理想椭圆上的点到两焦点的距离之和为常数(椭圆长轴2a). 实际中, 由于不可避免地存在误差, D与2a的偏差(即残差)为 2.RDa=-(3) 将上式应用于所有边缘点得到一系列的Ri {i = 1,2,,N}. 按最小二乘法, 要求 2 1 . N i i RMin = =(4) 迭代拟合过程中, 我们设置如下约束准则: 剔除残差绝对值|Ri|2.6σ (σ 为最小二乘拟合后的单位权误差)的边缘点; 如相邻两次迭代计算的任一焦点坐标的绝对偏差均小于0.001像素则迭代结束. 实际中, 我们发现迭代过程通常只需3~5次即可收敛. 处理后的边缘见图2(d). 图2 土星光环边缘检测的图像处理 (a) 包围光环边缘的矩形区域; (b) 经四邻点模板检测后的光环边缘; (c) 改进后的光环边缘(掏空内部边缘点); (d)经迭代椭圆拟合后的边缘点; (e) 用双线性内插做中心对称的去晕处理后; (f) 增加一常数背景水准后改进了视觉效果 2.2 晕的处理 由于土星光环具有极好的对称性, 我们可以依对称性对晕进行处理. 这里, 我们对已经测量了土星光环几何中心的原始图像按测定的中心点作为对称中心实施对称相减. 由于孀¹称中心不一定为整数值, 导致每一整像素点关于对称中心的对称点通常为非整像素位置点. 为了得到非整像素位置点的灰度值, 我们对该点的4个邻点应用双线性内插计算其灰度值. 具体计算式可以参考任何一本图像处理教材(例如Castleman的著作[15]). 经上述晕的处理后背景灰度通常很低(零附近). 对于暗卫星位置的进一步测量这是不合适的. 因为去晕后区域的星象与背景将形成过高的对比度, 不利于卫星位置测量的人机交互处理(见图2(e)). 根据我们从前的经验, 采用修正矩方法可以获得好的卫星位置测量精度. 这一方法不受卫星所在区域平均背景水准高低的影响(参见文献[11]). 这提供给我们机会来减少上述对比度的差异. 具体地, 我们可以对对称处理后的所有像素增加一个常数灰度值. 这一常数灰度值可以在对称处理前由处理者自由选择某一合适区域来确定. 图2(f)为最后处理得到的图像. 从该图中, 我们看到去晕处理后的卫星, 特别是近主星土卫1(Mimas)土卫2(Enceladus)的对比度有了明显改善. 无疑新处理方法让我们获得更多的机会来测量暗卫星的位置. 3 CCD视场定标和位置测量结果 3.1 卫星的位置测量与CCD视场的定标 与文献[11]中的测量方法一样, 我们采用修正矩方法测量土卫2~6的位置. 而土卫1尽管图像处理后具有好得多的可见性, 但它毕竟很暗且非常靠近光环(再次参考图2(f)), 有时候,我们不得不用光标凭视觉手工测量其位置. 根据Shen 等人[16]的研究结果, 土星主要卫星的TASS1.7理论与当代其他卫星理论相比有着最好的内部符合, 他们建议单独采用其4颗主要卫星(土卫3~6)的理论位置来确定CCD视场的定标参数. 本文即采用土星的4颗主要卫星(土卫3~6)的测量位置与理论位置进行比较来测定CCD视场的比例尺和指向(具体方法可以参考文献[11]), 其中土星的理论位置采用JPLDE4051)计算. 卫星和土星的理论位置计算到地心位置, 并考虑到大气折射周日视差以及切平面中心投影的影响. 表2给出了3个晚上CCD视场的比例尺和指向的测定结果. 表2 CCD视场的比例尺和指向的测定结果 资料集比例尺(/像素)指向/() 0201010.373888170.8914 020103a0.373871173.4648 020103b0.373895176.3346 0201040.373863179.4788 表2中每个晚上因为有其他观测任务而调整了CCD视场的指向, 其中1月3日晚上, 因为机械故障, 观测过程中调整了CCD底座, 指向也发生了变化. 3.2 位置测量结果 表3给出了土卫1~5(Mimas, Enceladus, Tethys, Dione, Rhea)以及土星相对于土卫6 (Titan) 位置测量的结果. 由表3可见, 所有卫星位置资料中, 土卫3~5相对于土卫6有最好的测量精度. 土卫1和2测量的结果精度较低, 特别是土卫1. 但我们看到这里土卫2的结果精度接近了土卫3~5的精度. 此外, 土卫1可测量位置的资料数虽然最少, 但是与我们过去测量的资料数或其他观测者测量的资料数相比其可测资料的比例已经很高了. 表4给出了此次观测中土卫1相对于土卫3~6观测资料之和的百分比, 也列出了文献[11]及其他观测者发表的资料百分比. 结果明显表明, 新的图像处理方法可以让我们获得比以前多得多的测量土卫1位置的机会. 最令人高兴的是从表3中我们看到测得的土星位置与理论位置之间有着很好的符合(3个晚上所有资料平均O-C的绝对偏差不大于0.01), 而且测量精度与土卫3~5一样, 达到了目前国际上土卫位置测量的最好精度(0.04~0.06). 020103a和020103b两资料子集有相对较大标准误差的原因主要是1月3日晚有较差的宁静度. 3.3 土卫测量结果与其他观测结果的比较 表5给出了本次试验观测的结果与同类CCD观测结果的比较. 由表5可见, 此次观测的 表3 土星卫星相对于土卫6位置测量的结果(单位: mas) 平均(O-C) 资料子集卫星或土星资料数 赤经方向标准误差赤纬方向标准误差 020101S1-Mimas19-246189106189 S2-Enceladus27-2431-1448 S3-Tethys27-8321033 S4-Dione271427-643 S5-Rhea271415531 Saturn277313337 020103aS1-Mimas12-22193345208 S2-Enceladus213959270 S3-Tethys21-639-1552 S4-Dione21544-463 S5-Rhea21743463 Saturn211435056 020103bS1-Mimas8-29896-11164 S2-Enceladus251866-774 S3-Tethys25-2063-867 S4-Dione251255-1253 S5-Rhea25160071 Saturn25053-852 020104S1-Mimas53-317150-16182 S2-Enceladus544342-3732 S3-Tethys541128032 S4-Dione54130-1431 S5-Rhea54-227632 Saturn54-1924634 TotalS1-Mimas92-28815557223 S2-Enceladus1272355-2055 S3-Tethys127-242-245 S4-Dione127739-1046 S5-Rhea127338445 Saturn127-438846 表4 土卫1相对于土卫3~6观测资料之和的百分比 资料比例(百分比) Harper等1997[8]9.2 Qiao等1999[9]2.8 Vienne等2002[10]5.3 Peng等2002[11]8.3 本工作18.1 表5 与其他观测结果的比较 (基于土卫3~5相对于土卫6的结果) (单位: mas) 平均(O-C) 观测者资料数 赤经方向标准误差赤纬方向标准误差 Harper等1997[8]891-770190 Qiao等1999[9]*3818090 Vienne等2001[10]4077174262 Peng等2002[11]451140-240 本工作381340-345 * 该观测资料采用极坐标系统 结果与文献[11]的结果几乎有一致的精度. 赤纬方向有稍大误差的原因除4.2中所述的大气宁静度因素外, 短时间(1.5~2.5 s)的露光也可能是导致测量精度稍有下降的原因. 4 结论 本文成功试验出一种基于土星光环边缘检测的图像处理技术来精确测定土星位置的方法,并进一步实施了土星本体和光环产生的晕的处理. 利用中国科学院云南天文台1 m望远镜上观测获得的3个晚上共计127帧CCD图像的实测处理表明: 土星光环几何中心位置的测量具有与土星主要卫星(土卫3~6)位置测量同样的精度, 单次位置测量的精度优于0.05. 这一精度达到了目前国际上土卫观测的最好精度. 此外, 晕的处理让我们获得了比前人多得多的机会来测量离土星本体和光环很近的土卫1和土卫2的位置. 我们相信这种新的图像处理技术将对CCD子午环观测土星位置以及长焦距望远镜观测土星暗卫星的位置均有重要的推动作用. 此外, 类似的处理方法将对其他行星(如天王星, 海王星)及其卫星位置的天体测量具有借鉴意义. 致谢 感谢南京大学天文系张承志教授和中国科学院陕西天文台沈凯先研究员阅读本文初稿并提出宝贵建议. 参 考 文 献 1 Sinclair A T. A theory of the motion of iapetus. MNRAS, 1974, 169: 591~605 2 Sinclair A T. The orbits of tethys, dione, rhea, titan and iapetus. MNRAS, 1977, 180: 447~495 3 Pascu D, Schmidt R E. Photographic positional observations of saturn. Astron J, 1990, 99: 1974~1984 4 Stone R C, Harris F H. CCD positions determined in the international celestial reference frame for the outer planets and many of their satellites in 1995~1999. Astron J, 2000, 119: 1985~1998 5 Stone R C. Positions for the outer planets and many of their satellites, IV. FASTT observations taken in 1999~2000. Astron J 2000, 120: 2124~2130 6 Stone R C. Positions for the outer planets and many of their satellites, V. FASTT observations taken in 2000~2001. Astron J, 2001, 122: 2723~2733 7 Rapaport M, Teixeira R, Campion J F L, et al. Astrometry of pluto and saturn with the CCD meridian instruments of Bordeaux and valinhos. Astron Astrophys, 2002, 383: 1054~1061 8 Harper D, Murray C D, Beurle K, et al. CCD astrometry of saturn's satellites 1990~1994. Astron Astrophys Suppl Ser, 1997, 121: 65~69 9 Qiao R C, Shen K X, Liu J R, et al. 1994~1996 CCD astrometric observations of saturn's satellites and comparison with theories. Astron Astrophys Suppl Ser, 1999, 137: 1~5 10 Vienne A, Thuillot W, Veiga C H, et al. Saturnian satellites observations made in brazil during the 1995 opposition with an astrometric analysis. Astron Astrophys, 2001, 380: 727~733 11 Peng Q Y, Vienne A, Shen K X. Positional measuring procedure and CCD observations for Saturnian satellites. Astron Astrophys, 2002, 383: 296~301 12 Pascu D, Seidelmann P K, Schmidt R E, et al. Astrometric CCD observations of Miranda: 1981~1985. Astron J, 1987, 93: 963~967 13 Vienne A, Duriez L. TASS1.6: Ephemerides of the major saturnian satellites. Astron Astrophys, 1995, 297: 588~605 14 Duriez L, Vienne A. Theory of motion and ephemerides of Hyperion. Astron Astrophys, 1997, 324: 366~380 15 Castleman K R. Digital Image Processing. 2nd ed. NJ: Prentice Hall, 1996, 117~119 16 Shen K X, Dourneau G, Qiao R C, et al. An analysis on satellite calibration method for CCD astrometry of saturn's satellites.