傅夏青，彭青玉，馮　杰

(1. 暨南大學(xué) 計(jì)算機(jī)科學(xué)系，廣東 廣州　510632；2. 暨南大學(xué) 中法天體測(cè)量、動(dòng)力學(xué)與空間科學(xué)聯(lián)合實(shí)驗(yàn)室，廣東 廣州　510632)

?

基于幾何扭曲模型的CCD圖像扭曲校正殘差及分析*1

傅夏青1 ,2，彭青玉1 ,2，馮杰1,2

(1. 暨南大學(xué) 計(jì)算機(jī)科學(xué)系，廣東 廣州510632；2. 暨南大學(xué) 中法天體測(cè)量、動(dòng)力學(xué)與空間科學(xué)聯(lián)合實(shí)驗(yàn)室，廣東 廣州510632)

摘要：在已知圖像幾何扭曲模型的基礎(chǔ)上，對(duì)2013年2月5日使用云南天文臺(tái)2.4 m望遠(yuǎn)鏡觀測(cè)近地小行星(99942)Apophis資料中的58幅CCD圖像進(jìn)行了幾何扭曲校正，并對(duì)校正后的圖像測(cè)出的星像位置進(jìn)行了幾何扭曲殘差分析。結(jié)果發(fā)現(xiàn)，星像位置X、Y兩個(gè)方向的幾何扭曲殘差約在 ± 0.02 pixel內(nèi)。為了探究影響幾何扭曲殘差大小的因素，設(shè)計(jì)了仿真實(shí)驗(yàn)。將模擬星像函數(shù)經(jīng)過設(shè)定的扭曲函數(shù)模型添加幾何扭曲，得到受扭曲的星像函數(shù)，并對(duì)扭曲前后的星像進(jìn)行幾何扭曲殘差分析。結(jié)果顯示，星像的半峰全寬(FWHM)越大，幾何扭曲模型的扭曲變化越快，幾何扭曲殘差越大。

關(guān)鍵詞：天體測(cè)量；幾何扭曲；扭曲校正；仿真

在天文觀測(cè)中，望遠(yuǎn)鏡光學(xué)成像系統(tǒng)經(jīng)常存在幾何扭曲，導(dǎo)致觀測(cè)得到的天體位置存在幾何扭曲誤差。如哈勃空間望遠(yuǎn)鏡(Habble Space Telescope, HST)的第二代廣域和行星照相機(jī)(Wide-Field Planetary Camera, WFPC2)拍攝的CCD圖像存在嚴(yán)重的幾何扭曲。其非常小視場(chǎng)(WFPC2視場(chǎng)中的寬場(chǎng)僅為80 × 80 arcsec2)的幾何扭曲最大可達(dá)5 pixel。這種由光學(xué)成像系統(tǒng)導(dǎo)致的星像位置移動(dòng)，就是幾何扭曲(Geometric Distortion, GD)。文[1]通過提高線性項(xiàng)的精度提出了一種改正幾何扭曲的方案，求解了哈勃空間望遠(yuǎn)鏡的幾何扭曲，使得哈勃空間望遠(yuǎn)鏡拍攝的圖像在行星相機(jī)視場(chǎng)精確到0.02 pixel，在寬場(chǎng)視場(chǎng)精確到0.01 pixel。解決了幾何扭曲問題后，哈勃空間望遠(yuǎn)鏡才真正發(fā)揮了其作為空間設(shè)備在天體測(cè)量上的潛能。在地面望遠(yuǎn)鏡CCD成像觀測(cè)中,也同樣存在視場(chǎng)幾何扭曲的影響。文[2-6]對(duì)地基天文望遠(yuǎn)鏡幾何扭曲的求解做了系統(tǒng)的研究，求解了云南天文臺(tái)2.4 m和1 m望遠(yuǎn)鏡不同濾光片的幾何扭曲模型并做了實(shí)際應(yīng)用。本文基于已知的幾何扭曲模型改正觀測(cè)的CCD圖像，并對(duì)其中產(chǎn)生的星像位置幾何扭曲殘差進(jìn)行了仿真研究。

1觀測(cè)資料和幾何扭曲模型

1.1觀測(cè)資料

本文所采用的圖像是2013年2月5日晚使用云南天文臺(tái)2.4 m望遠(yuǎn)鏡配備的云南暗弱天體光譜及成像儀(Yunnan Faint-Object Spectrog-raph and Camera, YFOSC)B濾光片觀測(cè)近地小行Apophis所得圖像中的58幀。因?yàn)樵颇习等跆祗w光譜及成像儀圖像的四周有黑框，所以先對(duì)原始圖像進(jìn)行適當(dāng)?shù)牟眉簦眉艉蟮膱D像大小為1 900 × 1 900像素。

1.2幾何扭曲模型

文[2-3]提出了一種簡(jiǎn)便有效的測(cè)定望遠(yuǎn)鏡CCD視場(chǎng)幾何扭曲的新方法，并應(yīng)用到云南天文臺(tái)1 m和2.4 m望遠(yuǎn)鏡抖動(dòng)觀測(cè)疏散星團(tuán)所得的圖像，求解了圖像的幾何扭曲模型。經(jīng)扭曲改正后，視場(chǎng)中星像的定位精度有了明顯提高[5-6]。本文使用的幾何扭曲模型正是用這種方法求解得到的。

如圖1，將大小為1 900 × 1 900像素的圖像劃分為19 × 19個(gè)大小為100 × 100像素的子區(qū)域，黑點(diǎn)表示子區(qū)域的中心像素點(diǎn)(格點(diǎn))。圖2為求解得到的云南天文臺(tái)2.4 m望遠(yuǎn)鏡B濾光片的幾何扭曲矩陣的部分?jǐn)?shù)值。矩陣大小為19 × 19，其中每一個(gè)數(shù)值代表對(duì)應(yīng)格點(diǎn)位置的幾何扭曲量。對(duì)于該幾何扭曲模型數(shù)值對(duì)應(yīng)的精度，文[2]做了相關(guān)研究，結(jié)果顯示絕大部分幾何扭曲量數(shù)值對(duì)應(yīng)的標(biāo)準(zhǔn)差小于0.05 pixel，并且各個(gè)標(biāo)準(zhǔn)差的分布是相對(duì)隨機(jī)的，與圖像坐標(biāo)沒有顯著的相關(guān)性。圖3是將幾何扭曲矩陣放大200倍后用矢量方式繪出。從圖3可以看到，幾何扭曲在圖像的上面角落位置比較大。

圖1格點(diǎn)

Fig.1Grids

圖2　幾何扭曲模型的數(shù)值矩陣形式

圖3　幾何扭曲模型的矢量圖

2幾何扭曲校正

(1)

由已知的幾何扭曲矩陣，圖像中任意位置處的幾何扭曲可由離該位置最近的4個(gè)格點(diǎn)的幾何扭曲量用雙線性插值的方法求出。

圖4幾何扭曲

Fig.4Geometric distortion

對(duì)圖像進(jìn)行扭曲校正的過程如圖5。黑點(diǎn)表示一幅圖像中不受扭曲的整像素點(diǎn)，受到幾何扭曲后偏移到了空心點(diǎn)的位置。如整像素點(diǎn)P(i,j)受到扭曲后移動(dòng)到P′(x,y)，那么用(2)式求出P′的位置。在扭曲圖像中插值求出P′位置的灰度值，并賦值給像素P。這樣就恢復(fù)了P像素點(diǎn)無扭曲的灰度值。對(duì)圖像中的每一個(gè)像素應(yīng)用此過程，就可以得到一幅扭曲校正后的圖像。

(2)

3實(shí)驗(yàn)結(jié)果

圖5像素扭曲示意圖

Fig.5Distorted pixel

對(duì)觀測(cè)的CCD圖像搜星得到存在幾何扭曲的星像坐標(biāo)，記為B(xgd,ygd)。用(3)式計(jì)算其扣除幾何扭曲后的星像坐標(biāo)C(x,y)。式中C的位置未知，其扭曲量也未知，所以無法直接計(jì)算C的位置。實(shí)際上相鄰位置的幾何扭曲幾乎相等(實(shí)際求解的幾何扭曲模型最大扭曲量約為±2 pixel)。所以先用B位置的幾何扭曲近似C位置的幾何扭曲，求出C的近似位置，然后再用迭代的方法計(jì)算C更精確的位置。用(4)式將兩組坐標(biāo)相減，得到星像的幾何扭曲。圖6中灰色點(diǎn)為星像位置的幾何扭曲分別隨星像位置的x、y坐標(biāo)的分布情況。

對(duì)扭曲校正后的CCD圖像搜星得到的星像中心，記為A(xgdc,ygdc)，與C(x,y)相減，得到扭曲校正后星像位置的幾何扭曲殘差，即(5)式。圖6中黑色點(diǎn)為星像位置的幾何扭曲殘差分別隨x、y坐標(biāo)的分布情況。將圖6中黑色點(diǎn)表示的幾何扭曲殘差量程縮小到± 0.05 pixel ，得到圖7。圖8為星像的幾何扭曲殘差放大200倍后的矢量分布圖。

(5)

圖6星像幾何扭曲和幾何扭曲殘差隨X/Y的變化

Fig.6The GD errors and residuals inX/Y

圖7幾何扭曲殘差隨X/Y的變化

Fig.7The GD residuals inX/Y

從圖6可以看到，圖像做扭曲改正前，星像的幾何扭曲最大可達(dá)1.7 pixel。扭曲改正后，星像的幾何扭曲殘差基本分布在0附近。從圖8可以看到，星像的幾何扭曲殘差分布的均值為0.008 pixel，標(biāo)準(zhǔn)差為0.009 pixel，并且與圖3相比明顯消除了幾何扭曲的區(qū)域系統(tǒng)性。在圖7中將殘差量程縮小后，可以看到幾何扭曲殘差大部分在± 0.02 pixel內(nèi)，但發(fā)現(xiàn)幾何扭曲殘差的分布呈現(xiàn)一定的系統(tǒng)趨勢(shì)。出現(xiàn)這個(gè)現(xiàn)象，可能的一個(gè)原因是幾何扭曲模型在圖像邊緣處的變化比較大，造成圖像改正后的幾何扭曲殘差也比較大。為了進(jìn)一步研究影響幾何扭曲殘差的因素，本文做了星像幾何扭曲的仿真實(shí)驗(yàn)。

4仿真實(shí)驗(yàn)

4.1仿真實(shí)驗(yàn)設(shè)計(jì)

圖8幾何扭曲殘差矢量圖

Fig.8Vectorgraph of GD residuals

如圖9，(a)為用一維連續(xù)高斯函數(shù)模擬的星像f(x)，表達(dá)式為(6)式，xc是模擬星像的中心。(b)為用正弦函數(shù)模擬的扭曲函數(shù)模型g(x)，表達(dá)式為(7)式。圖中為了突出星像扭曲的效果，在對(duì)扭曲函數(shù)繪圖時(shí)縮小了扭曲函數(shù)的周期。實(shí)際上，一個(gè)星像范圍內(nèi)的扭曲變化非常緩慢。(c)為對(duì)模擬星像函數(shù)添加泊松噪聲，再根據(jù)扭曲函數(shù)添加扭曲后得到的連續(xù)星像函數(shù)。根據(jù)扭曲前后模擬星像函數(shù)的積分面積不變，由(8)式可以推導(dǎo)出受扭曲后的連續(xù)星像函數(shù)h(t)，如(9)式。

根據(jù)文[7]對(duì)不同星像定心算法的比較和精度分析表明，高斯擬合法是一種精度相對(duì)較高的定心算法，且星像分布與高斯函數(shù)相似程度越高，則測(cè)量精度越高。實(shí)際求解的幾何扭曲模型最大扭曲量約為2 pixel，扭曲后的星像分布函數(shù)非常接近高斯函數(shù)。因此采用高斯擬合法對(duì)受扭曲后的星像定心。具體地，對(duì)受扭曲后的連續(xù)星像函數(shù)在其± 2.5σ像素范圍內(nèi)采樣(σ是高斯函數(shù)的標(biāo)準(zhǔn)差)。設(shè)此范圍內(nèi)覆蓋的像素區(qū)間為(i,i+n)。即對(duì)(i,i+1)、(i+1,i+2)…(i+n-1,i+n)區(qū)間分別積分，結(jié)果得到一組離散數(shù)據(jù)點(diǎn)，用圖9(c)中的點(diǎn)表示。對(duì)這些離散數(shù)據(jù)點(diǎn)用(10)式做最小二乘擬合，求出受扭曲后的星像中心xc′。實(shí)驗(yàn)中對(duì)每個(gè)像素區(qū)間積分時(shí)，采用兩點(diǎn)高斯數(shù)值積分[8]。

(10)

(6)～(10)式中，B為背景值，實(shí)驗(yàn)中設(shè)定星像扭曲前后背景值不變；H為星像亮度值；A是扭曲函數(shù)幅度，實(shí)驗(yàn)中參考真實(shí)資料求解得到的扭曲矩陣，將A設(shè)為2 pixel；T是扭曲函數(shù)的周期，T越小，扭曲變化越快；H′、xc′、σ′ 為對(duì)離散點(diǎn)做最小二乘擬合后求得的高斯函數(shù)的參數(shù)。

圖9(a)模擬的星像函數(shù); (b)扭曲函數(shù); (c)扭曲后的星像函數(shù)

Fig.9(a) Simulated star profile function; (b) GD function; (c) Distorted star profile function

4.2仿真實(shí)驗(yàn)結(jié)果

設(shè)圖像的寬度是1 024 pixel，在[0, 1023]區(qū)間內(nèi)隨機(jī)生成500個(gè)一維模擬星像函數(shù)。由(11)式計(jì)算幾何扭曲殘差，實(shí)驗(yàn)結(jié)果如圖9。

(11)

圖10中(a)表示設(shè)定扭曲函數(shù)的周期T=512 pixel，模擬星像的全峰半寬分別為13、9、4 pixel時(shí)，幾何扭曲殘差隨星像x坐標(biāo)的分布情況。表1是它的數(shù)值統(tǒng)計(jì)。(b)表示設(shè)定模擬星像的FWHM=9 pixel，扭曲函數(shù)的周期分別為256、512、1 024 pixel時(shí)，幾何扭曲殘差隨x坐標(biāo)的分布情況。表2是它的數(shù)值統(tǒng)計(jì)。

從圖10(a)和表1可以看到，當(dāng)扭曲函數(shù)的周期相同時(shí)，星像的全峰半寬越大，幾何扭曲殘差越大。例如在固定周期T=512 pixel情況下，當(dāng)FWHM=4 pixel時(shí)，幾何扭曲殘差最大為0.006 pixel，標(biāo)準(zhǔn)差為0.001 86 pixel；當(dāng)FWHM=9 pixel時(shí)，幾何扭曲殘差最大增至0.01 pixel，標(biāo)準(zhǔn)差增大到0.002 76 pixel；當(dāng)FWHM=13 pixel時(shí)，幾何扭曲殘差最大值雖然幾乎沒有增大，但標(biāo)準(zhǔn)差達(dá)到了0.004 pixel。從圖10(b)和表2可以看到，當(dāng)星像的全峰半寬相同時(shí)，扭曲函數(shù)的周期越小，幾何扭曲殘差越大。例如在固定FWHM=9 pixel情況下，當(dāng)T=1 024 pixel時(shí)，幾何扭曲殘差最大為0.007 5 pixel；當(dāng)T=256 pixel時(shí)，幾何扭曲殘差最大達(dá)到0.015 8 pixel。也就是說，星像的全峰半寬越大，扭曲函數(shù)的周期越小，幾何扭曲殘差越大。因?yàn)樾窍竦娜灏雽捲酱?，扭曲函?shù)的周期越小，都將導(dǎo)致一個(gè)星像范圍內(nèi)所包含的扭曲變化越大。從圖10還可以看到，幾何扭曲殘差的變化趨勢(shì)跟扭曲函數(shù)的變化一致。圖7中幾何扭曲殘差的分布在圖像的四周比較大，中間部分則較小。對(duì)應(yīng)圖3的幾何扭曲模型，幾何扭曲在圖像的角落處變化較快，相當(dāng)于仿真實(shí)驗(yàn)中的扭曲函數(shù)的周期T較小的情況；幾何扭曲在圖像的中間部分變化較小，相當(dāng)于仿真實(shí)驗(yàn)中扭曲函數(shù)的周期T較大的情況。從而導(dǎo)致了圖7中分布在圖像四周的星像幾何扭曲殘差較大，分布在圖像中間的星像幾何扭曲殘差較小。

圖10(a)相同周期下不同全峰半寬的扭曲殘差分布; (b)相同全峰半寬下不同周期的扭曲殘差分布

Fig.10(a) The GD residuals in the same period but with different FWHM; (b) The GD residuals with the same FWHM but in different periods

表1相同周期不同全峰半寬下扭曲殘差數(shù)值統(tǒng)計(jì)

Table 1The GD residuals statistics in the sameperiod but with different FWHM

FWHM(T=512)/pixelMax/pixelSD/pixel40.00590.0018690.01010.00276130.01030.00400

表2相同全峰半寬不同周期下扭曲殘差數(shù)值統(tǒng)計(jì)

Table 2 The GD residuals statistics in the same FWHM but in different periods

T(FWHM=9)/pixelMax/pixelSD/pixel2560.01580.005645120.00870.0027210240.00750.00262

5結(jié)論

在已經(jīng)求解出幾何扭曲模型的情況下，對(duì)CCD圖像進(jìn)行幾何扭曲校正。校正后，星像位置在X、Y兩個(gè)方向的幾何扭曲殘差約在± 0.02 pixel內(nèi)。對(duì)星像位置的幾何扭曲殘差進(jìn)行仿真研究，分別測(cè)試了星像在相同的扭曲函數(shù)周期下不同全峰半寬，以及星像在相同的全峰半寬下不同扭曲函數(shù)周期的幾何扭曲殘差的大小和變化規(guī)律。結(jié)果發(fā)現(xiàn)星像的全峰半寬越大，扭曲函數(shù)周期越小即扭曲模型變化越快，星像位置的幾何扭曲殘差越大。這個(gè)結(jié)論解釋了CCD圖像經(jīng)幾何扭曲校正后分布在圖像四周的星像其幾何扭曲殘差較大，而分布在圖像中間的星像其幾何扭曲殘差較小的現(xiàn)象。

參考文獻(xiàn)：

[1]Anderson J, King I R. An improved distortion solution for the Hubble Space Telescope′s WFPC2[J]. The Publications of the Astronomical Society of the Pacific, 2003, 115(803): 113-131.

[2]Peng Q Y, Vienne A, Zhang Q F, et al. A convenient solution to geometric distortion and its application to Phoebe′s observations [J]. The Astronomical Journal, 2012, 144(6): 170-179.

[3]Zhang Qingfeng, Peng Qingyu, Zhu Zi. Preliminary results of solving the problem of geometric distortion for the 2.4m telescope at Yunnan Observatory[J]. Research in Astronomy and Astrophysics, 2012, 12(10): 1451-1456.

[4]彭青玉, 涂兵. 求解CCD圖像幾何扭曲的初步結(jié)果[J]. 中國(guó)科學(xué): 物理學(xué) 力學(xué) 天文學(xué), 2011, 41(9): 1126-1130.

Peng Qingyu, Tu Bing. Preliminary results of solving geometric distortions for a CCD image[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2011, 41(9): 1126-1130.

[5]楊純一, 彭青玉, 張慶豐. 2.4米望遠(yuǎn)鏡天體測(cè)量試驗(yàn)[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2013, 10(2): 156-161.

Yang Chunyi, Peng Qingyu, Zhang Qingfeng. An astrometric experiments on the 2.4m telescope of the Yunnan Observatory[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2013, 10(2): 156-161.

[6]Peng Q Y, Wang N, Vienne A, et al. Precise CCD positions of Phoebe in 2011-2014[J]. Monthly Notices of the Royal Astronomical Society, 2015, 449(3): 2638-2642.

[7]李展, 彭青玉, 韓國(guó)強(qiáng). CCD圖像數(shù)字定心算法的比較[J]. 天文學(xué)報(bào), 2009, 50(3): 340-348.

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

[8]施妙根, 顧麗珍. 科學(xué)和工程計(jì)算基礎(chǔ)[M]. 北京: 清華大學(xué)出版社, 1999.

*基金項(xiàng)目：國(guó)家自然科學(xué)基金 (11273014, U1431227) 資助.

收稿日期：2015-09-39；

修訂日期：2015-11-23

作者簡(jiǎn)介：傅夏青，女，碩士. 研究方向：圖像處理. Email: cherishfu2015@sina.com

中圖分類號(hào)：P123.2+2

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào):1672-7673(2016)03-0382-07

Residuals and Analysis of Geometric Distortion Correction of CCD Images Based on Geometric Distortion Pattern

Fu Xiaqing1,2, Peng Qingyu1,2, Feng Jie1,2

(1. Department of Computer Science, Jinan University, Guangzhou 510632, China, Email: 897533236@qq.com 2. Sino-French Joint Laboratory for Astrometry, Dynamics and Space Science, Jinan University, Guangzhou 510632, China)

Abstract:Fifty eight CCD frames taken by the 2.4m telescope of Yunnan Observatories on Feb 5, 2013 were corrected for their Geometric Distortions (GDs) with the solved GD pattern by the approach proposed in Reference[2]. It is found that the GD residuals of the pixel position of stars still exist about ± 0.02 pixels in both X and Y directions while they are bigger at some corners and smaller in the middle part of each CCD frame. A simulation experiment was then carried out to analyze the cause that affects the GD residuals. Our results show that the GD residuals grow bigger when the FWHM (Full Width at Half Maximum) of a CCD stellar image becomes bigger or the GD pattern changes fast. This experiment can explain the variation of the GD residuals of stars after correcting the GDs.

Key words:Astrometry; Geometric distortion; Distortion correction; Simulation

CN 53-1189/PISSN 1672-7673