張 介，白金明，范玉峰

(1. 中國科學院云南天文臺，云南 昆明 650011; 2. 中國科學院研究生院，北京 100049;3. 中國科學院天體結構與演化重點實驗室，云南 昆明 650011)

云南天文臺麗江觀測站2.4 m望遠鏡是我國口徑最大的通用型光學望遠鏡，其主鏡直徑2.4 m，望遠鏡的主體由液壓支撐，主鏡由高壓氣體支撐，望遠鏡和圓頂?shù)尿寗雍透鞣N修正及報警等由計算機全自動控制，可支持遠程操作和自動操作。望遠鏡的成像質量好，分辨率、指向精度和跟蹤精度自動化程度都較高，綜合性能在國際上同級望遠鏡中處于中上水平。

2.4 m望遠鏡的開環(huán)跟蹤精度達到0.5″/1.4″/2″(10/30/60 min)，啟用自動導星系統(tǒng)時，其閉環(huán)跟蹤精度可以達到0.5″(60 min)，較好地保證了云南天文臺暗天體相機及分光儀(Yunnan Faint Object Spectrograph and Camera, YFOSC)、麗江地外行星探測儀(LiJiang Exoplanet Tracker, LiJET)等觀測終端進行長曝光時觀測數(shù)據(jù)質量。望遠鏡的自動導星系統(tǒng)采用偏置導星方式，在主光路內卡焦焦點前加入45°反射鏡提取40′×40′全視場內約4′×4′區(qū)域的視場作為導星CCD的視場，然后通過跟蹤導星區(qū)域內星像實現(xiàn)對主視場星像的跟蹤，導星鏡可以在與主光軸垂直的徑向調整位置，一般在最外側，不提取主視場10′×10′的視場，具有體積小、產熱量低、幾乎不影響鏡面視寧度的優(yōu)點。然而，目前的自動導星系統(tǒng)是制造商TTL為其早期制造的望遠鏡研制的，軟硬件都存在一些缺陷，運行不太穩(wěn)定，時常出現(xiàn)故障，中斷望遠鏡的自動跟蹤，影響了觀測數(shù)據(jù)的質量。參考同為TTL公司的2 m級英國利物浦望遠鏡自動導星系統(tǒng)升級經(jīng)驗，需要升級麗江2.4 m望遠鏡現(xiàn)有自動導星系統(tǒng)，以進一步提高現(xiàn)有望遠鏡系統(tǒng)的導星精度和效率。本文針對這一需求，開發(fā)系統(tǒng)中的自動導星定心模塊，用于后續(xù)的自動導星系統(tǒng)開發(fā)。

該自動導星系統(tǒng)首先依據(jù)望遠鏡觀測目標及導星系統(tǒng)的可觀測范圍，從GSC哈勃導星星表查詢可以用于導星的最佳星目標，接著通過移動導星反射鏡位置將最佳星目標引入導星CCD視場，然后連續(xù)地從自動導星CCD獲取導星圖像并計算星像位置偏移量，作為控制算法的誤差信號，指導望遠鏡的跟蹤修正。由于涉及較多的功能以及硬件接口，因此選用面向對象的C++語言開發(fā)軟件系統(tǒng)，并選用較成熟的開源軟件代碼庫，目標開發(fā)出一套有很強拓展性的相對獨立的自動導星系統(tǒng)。

本文主要介紹如何利用機器視覺方法從CCD圖像中提取用于導星的星像及計算星像中心實際位置。第1節(jié)介紹自動導星定心系統(tǒng)及通用的自動導星定心算法；第2節(jié)介紹自動導星系統(tǒng)星像提取算法及相關參數(shù)估計方法；第3節(jié)首先介紹利用LM算法非線性最小二乘擬合對星像二維高斯擬合，計算實際星像中心位置，然后使用2.4 m望遠鏡的觀測數(shù)據(jù)對系統(tǒng)進行測試，最后將高斯擬合結果與IRAF軟件二維高斯擬合結果相比較；第4節(jié)說明了星像偏移量的六常數(shù)模型計算方法。

1 自動導星定心系統(tǒng)及定心算法介紹

自動導星定心系統(tǒng)工作流程如圖1(a)，其軟件界面如圖1(b)，系統(tǒng)以Linux(Debian)[1]為平臺，利用開源代碼庫WXWIDGETS[2-3]開發(fā)圖形用戶界面(Graphical User Interface, GUI)，CFITSIO庫讀寫FITS文件(或SBIG相機讀寫驅動)，使用計算機視覺庫OPENCV/LAPACK[4-5]開發(fā)星像識別算法，MATHGL[6]庫生成各種圖像以及LEVMAR[7]庫作最小二乘擬合，算法底層為LAPACK矩陣計算，可以快速完成OPENCV和LEVMAR程序中所需矩陣運算，LAPACK有很強的拓展能力，可以實現(xiàn)多線程或者多臺計算機的并行計算，大幅提高運算效率。系統(tǒng)現(xiàn)可以從FITS(Pence[8])文件讀入數(shù)據(jù)、系統(tǒng)自動統(tǒng)計完成閾值設定、星像中心計算等。同時提供界面可以監(jiān)視CCD圖像，調試參數(shù)及相關結果統(tǒng)計顯示。

圖1 自動導星系統(tǒng)工作流程及自動導星軟件圖形用戶界面

Fig.1 The workflow of the auto-guiding system and the software Graphical User Interface of the system

一般的定心算法[9-14]主要有閾值一階矩質心、平方加權質心、高斯擬合中心[15-16]及點擴散函數(shù)相關運算質心算法[17]。若定義圖像中坐標(x,y)的修正灰度值為G(x,y)，則一階矩質心為：

類似地考察圖像灰度值的對稱性，可以計算平方加權質心：

如果一階矩質心與平方加權質心之間各個方向誤差隨機分布，說明CCD中天光背景均勻，反之需要天光背景補償。

高斯擬合公式：

高斯擬合算法的Stone[13,18]簡化公式：

高斯公式用于擬合的對數(shù)形式：

式中，B為背景天光值；P為星像最大灰度值；H對應實際星像峰值；R為擬合的高斯分布標準差。測試表明，完整模型比Stone簡化模型擬合有更好的擬合精度。

此外,系統(tǒng)中所用橢圓擬合算法，使用OPENCV庫函數(shù)效率極高，但返回結果為整型數(shù)據(jù)，有舍入誤差，只能用于確定星像區(qū)域，中心位置不可靠。

2 機器視覺星像搜索算法與參數(shù)估計

OPENCV提供了一系列機器視覺處理算法，這些算法基于LAPACK矩陣計算庫可以快速高效完成星像輪廓識別。針對CCD圖像中星像集中的特點，首先對8位圖像中值濾波去除異常噪聲，接著對圖像域值濾波并二值化，通過設定高于天光背景的域值可以區(qū)分星像與背景，然后可用Canny邊緣[19]檢測得到星像輪廓，最后再用最小二乘法擬合橢圓輪廓，估計星像區(qū)域。為簡化運算及數(shù)據(jù)存貯量，使用8位對數(shù)化整數(shù)數(shù)據(jù)完成星像識別，再使用32位原數(shù)據(jù)計算星像中心位置。

2.1 Canny邊緣檢測、圖像濾波與橢圓擬合

Canny邊緣檢測的算法是集低通濾波與邊界檢測于一體的算法，其內容如下：

第一步：利用高斯算子對圖像平滑卷積濾波去除噪聲，再計算圖像各點的灰度變化梯度，實際OPENCV的Canny算法源程序中直接采用同時有高斯平滑和邊緣檢測效果的Sobel算子與原圖像卷積計算，分別計算X與Y方向一階圖像差分。計算中采用3×3的模板計算。其數(shù)值如下：

第二步：計算絕對值范數(shù)或L2范數(shù)作為梯度強度和計算梯度方向。

梯度方向θ=arctan(Dy/Dx)，其中Dy與Dx是由Sobel算法算得的Y與X方向一階差分值。

第三步：梯度圖像非極大值抑制，將非局部最大梯度值點設為零。

第四步：雙閾值檢測和連接邊緣，沿梯度方向將圖像中梯度強度大于高閾值的都存為邊界點，低于高閾值且高于低閾值的梯度強度保留，再利用連通性篩選保留的梯度坐標，將與大于高閾值邊界點連通的梯度強度保留，其余無效區(qū)域設為零。

對圖像的濾波包括中值濾波去除異常數(shù)據(jù)點，但濾波僅對提取輪廓的8位數(shù)據(jù)進行，對用于擬合及定心計算的32位數(shù)據(jù)不進行濾波。圖像Canny提取輪廓圖像前利用天光背景估計值作為閾值將圖像二值化，可以大幅降低Canny算法對梯度強度閾值的敏感性，程序中使用(50, 125)的閾值對可以獲得比較好的效果。

然后利用輪廓提取函數(shù)可以獲得相互隔離又獨立連通的輪廓，最后使用最小二乘法擬合橢圓輪廓，獲得星像區(qū)域。圖2中測試文件為2.4 m望遠鏡觀測數(shù)據(jù)，視場10′×10′，曝光時間30 s，JR濾波片。

圖2 閾值154、160提取的星像輪廓及閾值155擬合的橢圓

Fig.2 The profiles of stellar images extracted with the threshold values 154 and 160,together with the ellipses fitted with the threshold value 155

2.2 圖像灰度值直方統(tǒng)計與CCD圖像天光背景估計

天光背景值在高斯擬合中具有很大影響，CCD圖像像素直方圖統(tǒng)計與實驗表明，直方圖中峰值對應灰度值為眾數(shù)，可以作為天光背景值。通?？梢詫?位對數(shù)化灰度直方統(tǒng)計峰值對應灰度值加2以上值判為有效星光灰度值。

程序中考慮到眾數(shù)附近的灰度值分布比較稠密，因此，在圖像的不同區(qū)域疊加后只有分布在眾數(shù)附近的灰度值才有可能在同子位置的值非常相近。程序首先將原數(shù)據(jù)做3×3小鄰域的均值濾波，獲得每個點鄰域內的均值作為該點的新灰度，并將濾波結果圖像邊緣近1/10 CCD尺寸的區(qū)域屏蔽，接著將有效區(qū)域分為11×11不重疊的子區(qū)域，然后對相鄰子區(qū)域的灰度值對應相減取絕對值，再與其他相鄰子區(qū)域計算結果重疊相加，接著求取最小值所在位置作為天光背景的眾數(shù)所在位置，最后將各個子區(qū)域此對應位置的灰度值相加取平均，作為天光背景值的眾數(shù)估計，同時也是天光背景值估計。同理地利用均值濾波對3×3鄰域標準差計算得到標準差，再求得區(qū)域疊加的最小位置，求得天光背景標準差σ(眾數(shù))估計，當然算法對星像過度稠密及天光背景不均勻的圖像可能會有估計不準的問題，目前算法還未對各星像分別計算天光背景值。

圖3 (a) 3×3鄰域均值的區(qū)域疊加圖像;(b) 3×3鄰域標準差的區(qū)域疊加圖像

Fig.3 (a) The image masked with a 3×3 mean-value filter;(b) The standards deviations of the masked image

為了在橢圓輪廓內限制有效區(qū)域，選擇灰度值接近1/5峰值處為邊界，對應于IRAF中測光孔徑值，可以減小天光背景的影響，獲得較好的擬合結果。

3 高斯函數(shù)擬合算法

正常星像受大氣影響，圖像灰度分布近似于二維高斯分布，且在各個方向應當有相同的標準差。系統(tǒng)采用非線性函數(shù)最小二乘擬合的方法，將參數(shù)估計問題轉化為最小化目標函數(shù)問題，得到獨立噪聲干擾下中心位置和標準差的最大似然估計。利用Levenberg-Marquardt優(yōu)化算法[9,20]擬合星像，該算法以均方誤差為目標函數(shù)，兼有梯度下降和牛頓-高斯方法的下降的速度，不直接求取復雜Hessian矩陣，用Jacobi行列式估計擬Hessian矩陣，程序中采用LEVMAR開源代碼作高斯最小二乘擬合。

為避免復雜的梯度函數(shù)，采用對數(shù)化數(shù)據(jù)擬合，全高斯公式Jacobi行列式為：

圖4 高斯擬合結果

Fig.4 Gaussian-fitting results

圖5是對2.4 m米望遠鏡YFOSC觀測數(shù)據(jù)高斯擬合星像位置與IRAF軟件高斯擬合結果比較，其最大誤差不超過0.08 Pixel。擬合得到高斯分布標準差為2.5，與當時記錄的視寧度為1.7″相符。

圖5 麗江2.4 m望遠鏡數(shù)據(jù)高斯擬合結果

Fig.5 The residual errors of the stellar positions from the Gaussian fitting of the stellar profiles observed with the Lijiang 2.4m telescope of the Yunnan Observatory

4 星像偏移量計算

自動導星算法對同一天區(qū)、相鄰曝光時間的兩幅圖像分別計算星像中心位置，并對相對應的星像計算位置偏移量，作為自動導星的誤差信號。

一般認為由于相鄰CCD圖像間可能存在平移、旋轉和縮放影響，可以用Stone[12,21-22]的算法，以六常數(shù)線性變換方程表示星像位置變換關系，如下，實際中i標識的數(shù)據(jù)量遠多于變換系數(shù)個數(shù)6。

編號97圖像：拍攝時間：2011-10-08T15∶02∶42.859，曝光時間30 s。

編號102圖像： 拍攝時間：2011-10-08T15∶07∶20.998，曝光時間50 s。

編號105圖像： 拍攝時間：2011-10-08T15∶10∶34.361，曝光時間50 s。

解六常數(shù)模型得：

abc97_102=[0.9999;2.3e-05;-0.9409]def97_102=[-8.3e-06;1.0001;2.3297]abc97_105=[1.0000;4.3e-05;-2.5640]def97_105=[2.3e-05;0.9999;3.5062]abc102_105=[1.0000;1.9e-05;-1.6231]def102_105=[3.1e-05;0.9999;1.1767]

(其中abc97_102表示97號圖與102號圖間X方向坐標變換系數(shù)[a102_97，b102_97，c102_97]，以此類推。)

圖6為六常數(shù)模型擬合的殘差，其標準化殘差都小于0.06 pixel，對其殘差進行Kolmogorov-Smirnov正態(tài)性檢驗，得到殘差數(shù)據(jù)都服從0.05顯著性水平的正態(tài)分布，可以認為六常數(shù)擬合已得到比較優(yōu)化的結果。

圖6 六常數(shù)模型97號～102號、97號～105號擬合殘差

Fig.6 The residual errors of the fits with the model of six coefficients solved for the No.97 and No.102 images versus those solved for No.97 and No.105 images

5 總 結

通過自主開發(fā)的自動導星系統(tǒng)對導星目標的圖像分析與算法實踐表明，對天光背景的眾數(shù)估計與直方統(tǒng)計的方法可以有效設定星像灰度閾值用于提取星像輪廓確定星像位置，利用基于LAPACK的LM算法可以快速高斯擬合，系統(tǒng)同時提供了完整的定心算法程序與參數(shù)估計方法，高斯擬合結果與IRAF軟件高斯擬合結果相比較不大于0.08 pixel(YFOSC終端0.282 ″/pixel)約0.022″，可以用六常數(shù)模型計算導星誤差量。而實際導星CCD像元為15 μm/pixel，對應天空角0.23 ″/pixel，可以得到更好的結果。

從算法運算效率來看，一般麗江2.4 m自動導星系統(tǒng)工作時導星CCD積分時間在1 s至25 s之間，而導星算法包括輸出控制在毫秒量級，也就是說實際制約自動導星系統(tǒng)采樣時間間隔的是CCD積分過程，而本算法系統(tǒng)雖然復雜，但即使在一般筆記本電腦(CPU T2400 1.83 GHz、內存2 G RAM、操作系統(tǒng)Debian6.02)，且LAPACK運算未開啟多線程并行計算時，運算時間可以滿足要求，另外采用導星計算與圖形用戶界面顯示多線程運行，降低圖形用戶界面監(jiān)控顯示的頻率，也可更進一步提高系統(tǒng)的響應速度。當然整套自動導星系統(tǒng)軟件涉及很多方面，工程具有一定規(guī)模，需在測試和使用中進一步完善。

致謝：感謝云南天文臺麗江觀測站2.4 m望遠鏡全體工作人員，感謝王夷博同學提供的觀測數(shù)據(jù)及天文數(shù)據(jù)處理方法的細致指導，感謝范緒亮同學提供的觀測FITS文件。

[1]Debian operation system main page[EB/OL]. [2012-04-19]. http://www.debian.org/.

[2]Wxwidgets Main Page[EB/OL]. [2012-04-19]. http://www.wxwidgets.org/.

[3]王強, Smart J, Hock K, et al. 使用wxWidgets進行跨平臺程序開發(fā), 2006.

[4]OpenCV Main Page [EB/OL]. [2012-04-19]. http://opencv.willowgarage.com/wiki/.

[5]LAPACK Main Page[EB/OL]. [2012-04-19]. http://www.netlib.org/lapack/.

[6]MathGL-library for scientific data visualization[EB/OL]. [2012-04-19]. http://mathgl.sourceforge.net/.

[7]LEVMAR Main Page [EB/OL]. [2012-04-19]. http://www.ics.forth.gr/～lourakis/levmar/.

[8]Pence W. CFITSIO Quick Start Guide[EB/OL]. [2012-04-19]. http://stderr.org/doc/libcfitsio-doc/quick/.

[9]王海涌, 費崢紅, 王新龍. 基于高斯分布的星像點精確模擬及質心計算[J]. 光學精密工程, 2009, 17(7): 1672-1677.

Wang Haiyong, Fei Zhenghong, Wang Xinlong. Precise simulation of star spots and centroid calculation based on Gaussian distribution[J]. Optics and Precision Engineering, 2009, 17(7): 1672-1677.

[10]潘波, 楊根慶, 劉勇. 星點質心定位算法最優(yōu)門限研究[J]. 光學精密工程, 2008, 16(9): 1787-1792.

Pan Bo, Yang Genqing, Liu Yong. Study on optimization threshold of centroid algorithm[J]. Optics and Precision Engineering, 2008, 16(9): 1787-1792.

[11]林潤芝, 楊學友, 鄒劍, 等. 面向大尺寸檢測CCD圖像中心提取精度的研究[J]. 傳感器與微系統(tǒng), 2010, 29(12): 51-53.

Lin Runzhi, Yang Xueyou, Zou Jian, et al. Study on the center extraction precision of image photographed by CCD for large scale inspection[J]. Transducer and Microsystem Technologies, 2010, 29(12): 51-53.

[12]李展, 彭青玉, 韓國強. CCD圖像數(shù)字定心算法的比較[J]. 天文學報, 2009, 50(3): 340-348.

Li Zhan, Peng Qingyu, Han Guoqiang. Comparison of digital centering algorithms based on CCD images[J]. Acta Astronomica Sinica, 2009, 50(3): 340-348.

[13]季凱帆, 王鋒. CCD圖像的二維修正矩定心方法[J]. 天文學報, 1996, 37(1): 85-90.

Ji Kaifan, Wang Feng. Two-dimensional modified moment centering algorithm in CCD images[J]. Acta Astronomica Sinica, 1996, 37(1): 85-90.

[14]季凱帆, 宋謙, 曹文達. CCD圖像的一維定心方法[J]. 云南天文臺臺刊, 1996(4): 69-74.

Ji Kaifan, Song Qian, Cao Wenda. The one-dimention centering algorithms of CCD image[J]. Publications of the Yunnan Observatory, 1996(4): 69-74.

[15]張志淵, 彭青玉. ePSF擬合法與Gaussian擬合法的比較[J]. 天文研究與技術——國家天文臺臺刊, 2010, 7(2): 132-139.

Zhang Zhiyuan, Peng Qingyu. Comparison of fitting ePSF and fitting gaussian-functions[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2010, 7(2): 132-139.

[16]王嫻, 蘇成悅, 王婷婷. 基于有效點擴散函數(shù)的高精度測光[J]. 天文研究與技術——國家天文臺臺刊, 2011, 8(1): 58-63.

Wang Xian, Su Chengyue, Wang Tingting. High-precision photometry based on an effective point spread function method [J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2011, 8(1): 58-63.

[17]李春艷, 謝華, 李懷鋒, 等. 高精度星敏感器星點光斑質心算法[J]. 光電工程, 2006, 33(2): 41-44.

Li Chunyan, Xie Hua, Li Huaifeng, et al. Centroiding algorithm for high-accuracy star tracker[J]. Opto-Electronic Engineering, 2006, 33(2): 41-44.

[18]Stetson P B. Photograhpic stellar photometry with the pds microdensitomter[J]. The Astronomical Journal, 1979, 84(7): 1056-1066.

[19]鄧林華, 柳光亁, 黃文娟, 等. 邊緣檢測算子在望遠鏡控制系統(tǒng)中的應用研究[J].微計算機信息, 2010(25): 57-58.

Deng Linhua, Liu Guangqian, Huang Wenjuan, et al. Applied research of edge detection operators in telescope control system[J]. Microcomputer Information, 2010(25): 57-58.

[20]Madsen K, Nielsen H B, Tingleff O. Methods for non-linear least squares problems[M]. Informatics and Mathematical Modelling Technical University of Denmark, 2004: 58.

[21]Stone Ronald C, Monet David G, Monet Alice K, et al. The flagstaff astrometric scanning transit telescope and star positions determined in the extragalactic reference frame[J]. The Astronomical Journal, 1996, 111(4): 1721-1743.

[22]齊朝祥, 于涌, 唐正宏. 施密特巡天星表的現(xiàn)狀和展望[J]. 天文學進展, 2008, 26(4): 349-359.

Qi Zhaoxiang, Yu Yong, Tang Zhenghong. Schmidt sky survey catalog: current status and perspectives[J]. Progress in Astronomy, 2008, 26(4): 349-359.