蔣夢源，尹東山，李海波，高玉平

天頂筒測量世界時的CCD星圖處理方法

蔣夢源1,2,3，尹東山1,2,3，李海波4，高玉平1,2,3

（1. 中國科學(xué)院大學(xué)，北京 100049；2. 中國科學(xué)院 國家授時中心，西安 710600；3. 中國科學(xué)院 時間頻率基準重點實驗室，西安 710600；4. 海軍91710部隊，和龍 133506）

為實現(xiàn)UT1自主測量，中國科學(xué)院國家授時中心建立了基于新型數(shù)字照相天頂筒的UT1測量網(wǎng)絡(luò)。該新型天頂筒采用CCD（charge coupled device）代替?zhèn)鹘y(tǒng)膠片作為感光元件，以提高觀測效率和測量精度。本文分析了天頂筒CCD星圖噪聲特性，發(fā)現(xiàn)以高斯噪聲為主，對比了多種濾波算法，發(fā)現(xiàn)高斯濾波算法對星圖的降噪效果更好。根據(jù)實拍星圖灰度分布的星像區(qū)域遠小于背景噪聲區(qū)域、星像灰度值高于背景噪聲灰度值等特征，選取合適閾值對整幅星圖進行分割，很好地將星像從背景噪聲中分離出來。比較了多種星像質(zhì)心位置的計算方法，發(fā)現(xiàn)加權(quán)質(zhì)心算法效果更好。通過實測數(shù)據(jù)檢驗，建立了一套適用于天頂筒星圖處理的方案。

CCD星圖;濾波;閾值分割;質(zhì)心提取

0 引言

天頂筒測量的基本工作原理是對天頂區(qū)域進行拍照，對拍攝到的星圖進行處理，提取星圖上星像位置，根據(jù)星圖上星像位置和恒星視位置，對星圖進行識別，得到星像與恒星的對應(yīng)關(guān)系，解算星圖中心點在天球上的指向點的視赤經(jīng)和視赤緯，即鉛垂線在天球上的指向點的視赤經(jīng)和視赤緯。

新型天頂筒采用電荷耦合器件（charge coupled device，CCD）進行成像，電荷耦合器件是一種新型半導(dǎo)體器件。它先將光學(xué)信號轉(zhuǎn)化為電信號，再根據(jù)模數(shù)轉(zhuǎn)化，將拍攝到的星圖以數(shù)字信號的形式進行存儲。相比于傳統(tǒng)底片，CCD存在畸變小、輕便、噪聲低、可集成化等優(yōu)點[1]。采用CCD作為感光元件來代替?zhèn)鹘y(tǒng)膠片，可以提高觀測效率和測量精度。

為實現(xiàn)UT1自主測量，中國科學(xué)院國家授時中心建立了基于新型數(shù)字照相天頂筒的UT1測量網(wǎng)絡(luò)。天頂筒拍攝的星圖以FITS（flexible image transport system）格式進行存儲，且像素為3 056×3 056，像元尺寸為9 μm×9 μm。FITS格式是一種專門用于天文工作的圖像存儲格式，它產(chǎn)生于20世紀70年代后期，它能很好地標準化，移植性好，并且可以包含關(guān)于圖像的描述參數(shù)，幾乎現(xiàn)在的天文軟件均支持這種格式[2]。

星圖處理是數(shù)據(jù)處理的關(guān)鍵部分。星圖處理提高圖像的信噪比，并提高最終的星像質(zhì)心位置測量精度。由于星圖存在噪聲，需要先對星圖進行濾波處理。并將星像所在區(qū)域與背景區(qū)域分割開來。通常一顆星像所占據(jù)的像素從幾個到上百個不等，而我們需要實現(xiàn)的是亞像素級別的位置測量精度，需要根據(jù)星像的灰度分布特征采用合適的算法進行計算。綜上，我們的星圖處理流程為：濾波、圖像分割和質(zhì)心位置計算。

1 CCD星圖處理方案

在天文圖像處理中常用的濾波算法有高斯濾波、均值濾波、中值濾波等，圖像分割算法有基于閾值分割算法、基于區(qū)域分割算法、基于邊緣分割算法等[3]，CCD相機自身結(jié)構(gòu)、觀測環(huán)境、大氣條件等因素會使得CCD星圖的噪聲類型、灰度直方圖和星像畸變等方面存在其自身特殊性。因此，對于CCD星圖處理，需要結(jié)合觀測臺站的實際情況，分析CCD星圖的特征，選擇適合CCD星圖的濾波方法和圖像分割方法。至于星像質(zhì)心位置計算，還需要根據(jù)我們提取到的星像信息具體分析，選擇位置誤差更小的方法。由于高性能CCD相機拍攝的星圖上星像較多，且每晚拍攝的底片有數(shù)百張之多，因此我們希望實現(xiàn)星圖的自動化處理，形成一套適用于實際星圖處理的方案。

1.1 濾波

濾波是圖像處理的第一步，目的是降低星圖噪聲。由于傳感器的材質(zhì)、觀測環(huán)境和元器件結(jié)構(gòu)等影響，CCD照相機拍攝到的星圖存在多種噪聲，如暗電流噪聲、電荷轉(zhuǎn)移噪聲、信號起伏噪聲、復(fù)位噪聲等[3-8]。根據(jù)噪聲分布特點，可以分為高斯噪聲、椒鹽噪聲和泊松噪聲。儀器電路各元器件自身噪聲及互相影響、圖像傳感器長期工作溫度過高等會在星圖上產(chǎn)生高斯噪聲。泊松噪聲即散粒噪聲，光子散粒噪聲符合泊松分布[9]。椒鹽噪聲又稱脈沖噪聲，它會隨機改變一些像素值[10]。

針對這幾類噪聲，常見的濾波方法有高斯濾波、均值濾波、中值濾波等。在圖像處理中，高斯濾波一般有兩種實現(xiàn)方式，一是用離散化窗口滑動卷積，另一種通過傅里葉變換。離散化窗口滑動卷積是選取一定尺寸大小的濾波器模板，依次遍歷圖像上的每個像素。高斯濾波就是將濾波器模板內(nèi)的像素灰度值做加權(quán)平均。均值濾波是對濾波器模板內(nèi)像素的灰度值做簡單平均計算，其缺點是不能很好地保護圖像細節(jié)，有可能會使圖像邊緣模糊化。中值濾波是采用濾波器模板內(nèi)像素灰度值的中位數(shù)代替濾波器模板中心點對應(yīng)的像素的灰度值，在平滑脈沖噪聲方面非常有效，同時它可以保護圖像尖銳的邊緣，是一種有效抑制椒鹽噪聲的濾波算法[10]。

為得到?jīng)]有星像的噪聲圖，使用洛南觀測站天頂筒在關(guān)閉圓頂、遮蓋鏡頭的情況下進行連續(xù)拍照，曝光時長0.5 s。以10幅圖像為一組，將一組圖像相同位置的灰度值求平均后得到平均后的圖像，分析其灰度分布直方圖。

圖1為平均后圖像的噪聲灰度分布直方圖，圖中曲線為高斯擬合曲線。由該圖中的直方圖及高斯擬合曲線的符合程度可看出，圖像噪聲的主要成分是高斯噪聲，這一類型的噪聲對應(yīng)的像素個數(shù)最多，灰度值也最大，對圖像的影響也最明顯。

圖1 灰度直方圖

對平均后的圖像分別計算其各列和各行灰度值的平均值，分析圖像灰度值的空間分布情況。

圖2和圖3分別表示平均后的圖像各行和各列的灰度平均值，由圖2可以看出各行灰度值的平均值在一個灰度值內(nèi)波動，由圖3可以看出列灰度值的平均值總體較穩(wěn)定，圖像兩側(cè)邊緣的列方向的灰度值的平均值較高，中間列中有個別列的灰度值平均值較小。各行各列的灰度值平均值總體穩(wěn)定，且不隨行列有明顯趨勢變化。

圖2 各行灰度值均值

圖3 各列灰度值均值

圖4和圖5分別表示平均后圖像各行和各列的灰度值的標準差，各行各列的灰度值的標準差總體平穩(wěn)，且不隨行列有明顯趨勢變化。

濾波器模板的大小一般是奇數(shù)，如3×3，5×5或者7×7，濾波器模板半徑越大，濾波耗時越長[11]，綜合考慮天頂筒所配置的計算能力等因素，決定采用3×3模板分別進行高斯濾波和中值濾波，并根據(jù)峰值信噪比指標對這兩種方法進行比較，選取合適的濾波方法。

圖4 各行灰度值標準差

圖5 各列灰度值標準差

1.2 圖像分割

圖像分割的目的是將星像區(qū)域從背景中分割出來。主要方法有：閾值分割法、基于區(qū)域的分割方法、基于邊緣的分割方法等[3]。

根據(jù)實拍星圖的像素灰度分布特性，星圖上星像和背景噪聲的灰度值處于不同灰度級范圍，星像的灰度值大于噪聲點的灰度值，且星像區(qū)域遠小于背景噪聲區(qū)域。

基于閾值的分割方法具有實現(xiàn)簡單、計算量小、性能較穩(wěn)定等特點[13]，且它適用于目標和背景處于不同灰度級范圍的圖像，利用圖像中需要提取的目標與背景在灰度上的差異，通過設(shè)置合理閾值實現(xiàn)目標與背景的分離。因此，選取基于閾值的分割方法。

閾值分割方法可表示為[14]

閾值分割中的閾值選取至關(guān)重要，如果閾值選取過小，會導(dǎo)致將很多噪聲誤認為是星像，如果閾值過大，也有可能剔除掉很多暗星。

閾值的選取有很多種方法，有最大類間方差法、迭代法等。每種方法都有其適用的圖像，根據(jù)圖像灰度特征選取不同的閾值選取法。

最大類間方差法是通過不斷調(diào)節(jié)閾值，使得通過此閾值分割的目標與背景之間的差異最大，進而認為此閾值選取最優(yōu)[15]。迭代法原理與最大類間方差法原理類似。

但是由于星像區(qū)域遠小于背景噪聲區(qū)域，利用最大類間法和迭代法計算出的閾值較小，閾值分割效果不好。當(dāng)目標與背景的面積大小比例懸殊時，閾值計算結(jié)果不準確，利用此閾值的圖像分割效果不好。

通常情況下，如果閾值分割得到的星像像素個數(shù)小于5，則認為該像素區(qū)域為噪聲，需要將其剔除。

1.3 質(zhì)心提取

星像質(zhì)心位置計算是為了從已經(jīng)提取到的星像區(qū)域上，得到高精度的星像質(zhì)心位置。

一個星像信號集中在近似圓形的區(qū)域內(nèi)，并且星像能量在區(qū)域內(nèi)呈近似高斯分布[16]。因此需要計算出星像區(qū)域內(nèi)近似能量最高的點，即星像質(zhì)心。

關(guān)于星像質(zhì)心提取的方法可分為基于灰度的質(zhì)心算法和基于邊緣的質(zhì)心算法[17]，其中基于灰度的方法有：質(zhì)心法、加權(quán)質(zhì)心法、高斯曲面擬合法等[18]。分別用這3種方法進行實驗，比對得出的星像質(zhì)心位置，選取最適合的質(zhì)心提取方法。

加權(quán)質(zhì)心法是質(zhì)心法的一種進化形式，加權(quán)質(zhì)心法和質(zhì)心法相比，其采用灰度值的平方代替了灰度值，它使得離中心較近且灰度值較大的像素點對于中心位置的計算貢獻更大。計算公式如式（6）和（7）所示：

因為星像在CCD上的成像為光斑，其灰度值近似服從高斯分布，因此可以用高斯曲面對其進行擬合。對二維高斯曲面公式兩邊求導(dǎo)，利用最小二乘法進行曲面擬合，為了計算簡便，也可以分別從和方向分別對星像進行擬合[19]。高斯曲線方程如下：

2 星像質(zhì)心提取結(jié)果

圖6為2018年11月23日23:01:56于洛南觀測站拍攝的一張原始星圖，曝光時長為0.5 s。使用上述分析的星圖處理方法，對星圖進行處理。

圖6 洛南站拍攝的星圖（2018-11-23）

圖7為高斯濾波后的星圖。高斯濾波后星圖的峰值信噪比為23.91，中值濾波后星圖的峰值信噪比為22.18，高斯濾波加中值濾波后星圖的峰值信噪比為22.80，高斯濾波后的星圖的峰值信噪比最大，說明高斯濾波算法更適合。

圖7 高斯濾波后星圖

圖8為閾值分割后的星圖。閾值分割得到的星像中，根據(jù)望遠鏡的拍攝能力，即使是最暗的星占據(jù)的像素個數(shù)也有6個以上，因此像素個數(shù)小于5的圖像看作噪聲點，需要將其剔除。

圖8 閾值分割后星圖

根據(jù)上述閾值分割提取到的星像，并將星像按照像素大小進行排序編號，序號越小的星像，其像素越小，序號越大的星像，其像素越大。利用3種方法分別進行星像質(zhì)心位置計算，并對3種方法的計算結(jié)果做對比。

圖9中點線表示高斯曲面擬合法和質(zhì)心法算出的質(zhì)心位置的行差，虛線表示高斯曲面擬合法和加權(quán)質(zhì)心法算出的質(zhì)心位置的行差，實線表示質(zhì)心法和加權(quán)質(zhì)心法算出的質(zhì)心位置的行差。圖10對應(yīng)圖9中的列差。

圖9 3種方法計算的星像質(zhì)心在行方向的差異

圖10 3種方法計算的星像質(zhì)心在列方向的差異

如圖9和10可以看出，實線條一直在0點附近，且波動較小。在橫坐標序號較小部分（對應(yīng)像素較小的星像），點線和虛線條在某些點波動較大。在橫坐標序號較大部分（對應(yīng)像素較大的星像），3線條一直在0點附近，且波動較小。無論星像像素大小，質(zhì)心法和加權(quán)質(zhì)心法得到的星像質(zhì)心位置相差很小。但是，高斯曲面擬合法得出的質(zhì)心位置和另兩種算法得出的質(zhì)心位置對于某些像素較小的星像差異較大。

為了分析高斯曲面擬合法和質(zhì)心法或者加權(quán)質(zhì)心法得到的結(jié)果在某些點上的差異較大，且在列方向上的差異更明顯的原因，我們挑選出這些在質(zhì)心列方向上3種方法計算結(jié)果差異較大的星像進行分析，發(fā)現(xiàn)這些星像區(qū)域較小，像素從幾個到十幾個不等，星像區(qū)域均為狹長區(qū)域，并且其灰度分布已不符合高斯分布。因此，在這些星像區(qū)域上，利用高斯曲面擬合法計算其質(zhì)心位置，誤差就會過大。與高斯曲面擬合法相比，加權(quán)質(zhì)心法的計算結(jié)果更穩(wěn)定，不存在誤差較大的情況。

由于天頂筒觀測不是跟蹤曝光，曝光期間的地球自轉(zhuǎn)會導(dǎo)致星像存在拖尾現(xiàn)象。除此之外，某些星像也存在畸變現(xiàn)象，畸變是由光學(xué)系統(tǒng)加工制作誤差導(dǎo)致的。

星圖上存在一定數(shù)目的星像區(qū)域狹長且灰度值分布不符合高斯分布的星像，這些星像并不適合用高斯曲面擬合法求其質(zhì)心位置。在星像提取過程中，我們已經(jīng)剔除了像素小于5的星像，如果我們再直接剔除掉這些星像，會導(dǎo)致星圖上星像過少，并對之后的計算帶來不便。因此我們需要一個適用于我們星圖上所有星像質(zhì)心位置提取，且計算結(jié)果穩(wěn)定的方法。

加權(quán)質(zhì)心法計算的星像質(zhì)心位置比較穩(wěn)定，且加權(quán)質(zhì)心法是質(zhì)心法的一種進化形式，其采用灰度值的平方代替了灰度值，它使得離中心較近且灰度值較大的像素點對于中心位置的計算貢獻更大。故采用加權(quán)質(zhì)心法計算星像質(zhì)心位置。

3 結(jié)論

[1] 劉美瑩. CCD天文觀測圖像的星圖識別和天文定位方法研究[D]. 西安: 中國科學(xué)院西安光學(xué)精密機械研究所, 2009.

[2] 冒蔚, 季凱帆, 李彬華, 等. CCD天體測量學(xué)[M]. 昆明: 云南科技出版社, 2003.

[3] 鄧廷權(quán), 焦穎穎. 圖像分割質(zhì)量評價的二型模糊集方法[J]. 計算機工程與應(yīng)用, 2011, 47(32): 217-220.

[4] 謝玉瑩. 星敏感器星圖仿真及星點提取算法評估研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 2017.

[5] 董雅瀟. 基于小視場星敏感器的星點提取與星圖識別算法研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2018.

[6] 劉美瑩. CCD天文觀測圖像的星圖識別和天文定位方法研究[D]. 西安: 中國科學(xué)院西安光學(xué)精密機械研究所, 2009.

[7] 蘇暢. 星敏感器星圖降噪預(yù)處理及其硬件實現(xiàn)[D].哈爾濱: 哈爾濱工業(yè)大學(xué), 2017.

[8] 李爽. 基于FPGA的CCD星敏感器星圖預(yù)處理[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 2007.

[9] 廖育富. 復(fù)雜工況背景下星圖預(yù)處理相關(guān)技術(shù)研究[D]. 成都: 中國科學(xué)院光電技術(shù)研究所, 2016.

[10] 周杰. 應(yīng)用于圖像處理的中值濾波改進算法[D]. 北京: 北京郵電大學(xué), 2007.

[11] 金龍, 王洪元, 張繼, 等. 實時DSP圖像處理高斯濾波優(yōu)化[J]. 制造業(yè)自動化, 2014, 36(12): 63-66.

[12] 王坤. 數(shù)字圖像分割和質(zhì)量評價方法的研究[D]. 沈陽: 東北大學(xué), 2006.

[13] 呂燕. 基于閾值算法圖像分割的研究[D]. 重慶: 重慶大學(xué), 2011.

[14] 聶僥. 基于CCD星敏感器的星圖識別[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2012.

[15] 蔡梅艷, 吳慶憲, 姜長生, 等. 改進Otsu法的目標圖像分割[J]. 電光與控制, 2007, 14(6): 118-151.

[16] 沈曉芳, 劉朝山. 星敏感器中的圖像預(yù)處理技術(shù)[J]. 河南科學(xué), 2016, 34(4): 471-476.

[17] SHORTIS M R, CLARKE T A, SHORT T. A comparison of some techniques for the subpixel location of discrete target images[J]. SPIE, 1994: 239-250.

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

[19] 張廣軍. 星圖識別[M]. 北京: 國防工業(yè)出版社, 2011.

CCD star map processing method for measuring UT1 with zenith telescope

JIANG Meng-yuan1,2,3, YIN Dong-shan1,2,3, LI Hai-bo4, GAO Yu-ping1,2,3

(1. University of Chinese Academy of Sciences, Beijing 100049, China; 2. National Time Service Center, Chinese Academy of Sciences, Xi’an 710600, China; 3. Key Laboratory of Time and Frequency Primary Standards, Chinese Academy of Sciences, Xi’an 710600, China; 4. Navy 91710, Helong 133506, China)

In order to realize the autonomous measurement of UT1, the National Time Service Center of Chinese Academy of Sciences established a Digital Photographic Zenith Telescope network. In order to improve the observing efficiency and measurement accuracy, the Digital Photographic Zenith Telescope uses CCD instead of the traditional film. In this paper, we analyzed the noise characteristics of the CCD star map, and found that Gaussian noise is the main noise source. So we compared a variety of filtering algorithms, and among them the Gaussian filtering algorithm exhibited better noise reduction performance than the others. According to the gray distribution characteristics of the real star map, the star image region is far smaller than the background noise region, and the star image gray value is higher than that of the background noise, we can select appropriate threshold value to segment the whole star map, which could separate the star image from the background noise well. In this paper, we compared several methods for calculating the centroid position of star image, and the weighted centroid algorithm was found to be more effective. Based on the research, we developed a set of optimal scheme for star map processing.

CCD (charge coupled device) star map; filter; threshold segmentation; weighted centroid algorithm

10.13875/j.issn.1674-0637.2020-02-0143-10

2019-10-23；

2020-01-21

陜西省自然科學(xué)基金資助項目（2018JM1031）；國家自然科學(xué)基金面上資助項目（11973046）；國家自然科學(xué)基金重大研究計劃資助項目（91736207）

蔣夢源，女，碩士，主要從事天體測量學(xué)研究。