导星系统论文：自动导星系统的改善设计探究范文

作者：张介白金明范玉峰单位：中国科学院云南天文台中国科学院大学中国科学院天体结构与演化重点实验室

本文主要介绍如何利用机器视觉方法从CCD图像中提取用于导星的星像及计算星像中心实际位置。第一节介绍自动导星定心系统及通用的自动导星定心算法，第二节介绍自动导星系统星像提取算法及相关参数估计方法。第三节首先介绍利用LM算法非线性最小二乘拟合对星像二维高斯拟合，计算实际星像中心位置，然后使用丽江2.4米望远镜的观测数据对系统进行测试，最后将高斯拟合结果与IRAF软件二维高斯拟合结果相比较。第四节说明了星像偏移量的六常数模型计算方法。

自动导星定心系统及定心算法介绍

自动导星定心系统工作流程如图1左图所示，其软件界面如图1右图所示，系统以Linux（Debian）(1)为平台，利用开源代码库WXWIDGETS(2,3)开发GUI界面、CFITSIO库读写FITS文件（或SBIG相机读写驱动）、使用计算机视觉库OPENCV/LAPACK(4)(5)开发星像识别算法、MATHGL(6)库生成各种图像以及LEVMAR(7)库作最小二乘拟合，算法底层为LAPACK矩阵计算，可以快速完成OPENCV和LEVMAR程序中所需矩阵运算，LAPACK有很强拓展能力，可以实现多线程或者多台计算机的并行计算，大幅提高运算效率。系统现可以从FITS(Pence(8))文件读入数据、系统自动统计完成阈值设定、星像中心计算等。同时提供界面可以监视CCD图像、调试参数及相关结果统计显示。一般的定心算法(10-14)主要有阈值一阶矩质心、平方加权质心、高斯拟合中心(15,16)，及(17)PSF相关运算质心算法。若定义图像中坐标(x,y)的修正灰度值为G(x,y),则:一阶矩质心为:如果一阶矩质心与平方加权质心之间各个方向误差随机分布，说明CCD中天光背景均匀，反之需要天光背景补偿。高斯拟合公式：高斯拟合算法的Stone(13,18)简化公式：及高斯公式用于拟合的对数形式：其中B为背景天光值，P为星像最大灰度值，H对应实际星像峰值，R为拟合的高斯分布标准差，测试表明，完整模型比Stone简化模型拟合有更好拟合精度。此外,系统中所用椭圆拟合算法，使用OPENCV库函数效率极高，但返回结果为整型数据，有舍入误差，只能用于确定星像区域，中心位置不可靠。

机器视觉星像搜索算法与参数估计

OPENCV提供了一系列机器视觉处理算法，这些算法基于LAPACK矩阵计算库可以快速高效完成星像轮廓识别。针对CCD图像中星像集中的特点，首先对8位图像中值滤波去除异常噪声，接着对图像域值滤波并二值化，通过设定高于天光背景的域值可以区分星像与背景，然后可用Canny边缘检测得到星像轮廓，最后再用最小二乘法拟合椭圆轮廓，估计星像区域。为简化运算及数据存贮量，使用8位对数化整数数据完成星像识别，再使用32位原数据计算星像中心位置。

1Canny边缘检测、图像滤波与椭圆拟合

Canny边缘检测(19)的算法是集低通滤波与边界检测于一体的算法，其内容如下：第一步：利用高斯算子对图像平滑卷积滤波去除噪声，再计算图像各点的灰度变化梯度，实际OPENCV的Canny算法源程序中直接采用同时有高斯平滑和边缘检测效果的Sobel算子与原图像卷积计算,分别计算X与Y方向一阶图像差分。计算中采用3x3的模板计算。其数值如下：第二步：计算绝对值范数或L2范数作为梯度强度和计算梯度方向。梯度方向arctan(Dy/Dx)，其中Dy与Dx是由Sobel算法算得的Y与X方向一阶差分值。

第三步：梯度图像非极大值抑制，将非局部最大梯度值点设为零。第四步：双阈值检测和连接边缘，沿梯度方向将图像中梯度强度大于高阈值的都存为边界点，低于高阈值且高于低阈值的梯度强度保留，再利用连通性筛选保留的梯度坐标，将与大于高阈值边界点连通的梯度强度保留，其余无效区域设为零。对图像的滤波包括中值滤波去除异常数据点，但滤波仅对提取轮廓的8位数据进行，对用于拟合及定心计算的32位数据不进行滤波。图像Canny提取轮廓图像前利用天光背景估计值作为阈值将图像二值化，可以大幅降低Canny算法对梯度强度阈值的敏感性，程序中使用（50，125）的阈值对可以获得比较好的效果。Fig.2Theoutlineofstarsimagesextractedwiththreshold154,160andtheellipseimagesfittedwiththreshold155然后利用轮廓提取函数可以获得相互隔离又独立连通的轮廓，最后使用最小二乘法拟合椭圆轮廓，获得星像区域。图2中测试文件为丽江2.4米望远镜观测数据，视场10′x10′,曝光时间30秒，JR滤波片，范绪亮同学提供。

2图像灰度值直方统计与CCD图像天光背景估计

天光背景值在高斯拟合中具有很大影响，CCD图像像素直方图统计与实验表明，直方图中峰值对应灰度值为众数，可以作为天光背景值。通常可以将8位对数化灰度直方统计峰值对应灰度值加2以上值判为有效星光灰度值。

程序中考虑到众数附近的灰度值分布比较稠密，因此，在图像的不同区域叠加后只有分布在众数附近的灰度值才有可能在同子位置的值非常相近。程序首先将原数据做3x3小邻域的均值滤波，获得每个点邻域内的均值作为该点的新灰度，并将滤波结果图像边缘近1/10CCD尺寸的区域屏蔽，接着将有效区域分为11x11不重叠的子区域，然后对相邻子区域的灰度值对应相减取绝对值，再与其他相邻子区域计算结果重叠相加，接着求取最小值所在位置作为天光背景的众数所在位置，最后将各个子区域此对应位置的灰度值相加取平均，作为天光背景值的众数估计，同时也是天光背景值估计。同理地利用均值滤波对3x3邻域标准差计算得到标准差，再求得区域叠加的最小位置，求得天光背景标准差σ（众数）估计，当然算法对星像过度稠密及天光背景不均匀的图像可能会有估计不准的问题，算法目前还未对各星像分别计算天光背景值。为了在椭圆轮廓内限制有效区域，选择灰度值接近1/5峰值处为边界，对应于IRAF中测光孔径值，可以减小天光背景的影响，获得较好的拟合结果。

高斯函数拟合算法

正常星像受大气影响，图像灰度分布近似于二维高斯分布，且在各个方向应当有相同的标准差。系统采用非线性函数最小二乘拟合的方法，将参数估计问题转化为最小化目标函数问题，得到独立噪声干扰下中心位置和标准差的最大似然估计。利用Levenberg–Marquardt优化算法(9,20)拟合星像，该算法以均方误差为目标函数，兼有梯度下降和牛顿－高斯方法的下降的速度，不直接求取复杂Hessian矩阵，用Jacobi行列式估计拟Hessian矩阵，程序中采用LEVMAR开源代码作高斯最小二乘拟合。

为避免复杂的梯度函数，采用对数化数据拟合，全高斯公式Jacobi行列式为：其中计算时P的初始值为星像峰值减天光背景值(+3σ)，(x,y)为图像中星像峰值坐标,为防止局部陷入给一定的初始偏差，测试表明算法有很好的收敛性。图4是高斯拟合结果与原数据比较，底部为残差密度图。下图是对丽江2.4米望远镜YFOSC观测数据高斯拟合星像位置与IRAF软件高斯拟合结果比较，其最大误差不超过0.08Pixel。拟合得到高斯分布标准差为2.5，与当时记录的视宁度为1.7角秒相符。

星像偏移量计算

自动导星算法对同一天区、相邻曝光时间的两幅图像分别计算星像中心位置，并对相对应的星像计算位置偏移量，作为自动导星的误差信号。一般认为由于相邻CCD图像间可能存在平移、旋转和缩放影响，可以用Stone(12,21,22)的算法，以六常数线性变换方程表示星像位置变换关系，如下，实际中i标识的数据量远多于变换系数个数6。以下是对同一目标的观测结果分析，其中：编号97图像：拍摄时间：2011-10-08T15:02:42.859，曝光时间30s。编号102图像：拍摄时间：2011-10-08T15:07:20.998，曝光时间50s。编号105图像：拍摄时间：2011-10-08T15:10:34.361，曝光时间50s。解六常数模型得：图6为六常数模型拟合的残差，其标准化残差都小于0.06pixel，对其残差进行Kolmogorov-Smirnov正态性检验，得到残差数据都服从0.05显著性水平的正态分布。可以认为六常数拟合已得到比较优化的结果。

总结

通过自主开发的自动导星系统对导星目标的图像分析与算法实践表明，对天光背景的众数估计与直方统计的方法可以有效设定星像灰度阈值用于提取星像轮廓确定星像位置，利用基于LAPACK的LM算法可以快速高斯拟合，系统同时提供了完整的定心算法程序与参数估计方法，高斯拟合结果与IRAF软件高斯拟合结果相比较不大于0.08pixel(YFOSC终端0.282″/pixel)约0.022角秒,可以用六常数模型计算导星误差量。而实际导星CCD像元为15µm/pixel,对应天空角0.23″/pixel,可以得到更好结果。

从算法运算效率来看，一般丽江2.4米自动导星系统工作时导星CCD积分时间在1s至25s之间，而导星算法包括输出控制在毫秒量级，也就是说实际制约自动导星系统采样时间间隔的是CCD积分过程，而本算法系统虽然复杂，但即使在一般笔记本电脑（CPUT24001.83Ghz、内存2GRAM、操作系统Debian6.02），且LAPACK运算未开启多线程并行计算时，运算时间可以满足要求，另外采用导星计算与GUI显示多线程运行，降低GUI监控显示的频率，也可更进一步提高系统响应速度。当然整套自动导星系统软件涉及很多方面，工程具有一定规模，需在测试和使用中进一步完善。致谢:感谢云南天文台丽江2.4米望远镜组全体工作人员在此期间的辛勤工作为文章完成提供了充裕时间，感谢王夷博同学提供的观测数据及天文数据处理方法的细致指导，感谢范绪亮同学提供的观测FITS文件。感谢WXWIDGETS(2,3)、OPENCV(4)、LAPACK(5)、MATHGL(6)、LEVMAR(7)开源软件作者的辛勤劳动。