曹西平,彭青玉,張慶豐
(1.暨南大學(xué)計(jì)算機(jī)科學(xué)系,廣東 廣州 510632;2.暨南大學(xué),中法天體測量、動(dòng)力學(xué)與空間科學(xué)聯(lián)合實(shí)驗(yàn)室,廣東 廣州 510632;3.廣東省高等學(xué)校光電信息與傳感技術(shù)重點(diǎn)實(shí)驗(yàn)室(暨南大學(xué)),廣東 廣州 510632)
根據(jù)Stone[1]的研究,在高斯函數(shù)擬合過程中為了避免像素相位誤差,星像在半高全寬(FWHM——Full Width at Half Maximum)下最好大于4 pix。然而,無論是地面成像還是空間成像,均可能存在星像在半高全寬下小于4 pix的情況。例如,在地面觀測中,美國海軍天文臺(tái)使用的雙筒照相儀得到的圖像就是欠采樣的,從他們發(fā)布的UCAC(The US Naval CCD Astrograph Catalog)系列星表中[2-4]可以看到。他們雖然采用了不同模型對(duì)CCD圖像中的星像進(jìn)行擬合[5],但像素相位誤差依然存在。在空間觀測中,Hubble空間望遠(yuǎn)鏡WFPC2照相機(jī)的欠采樣成像也是一個(gè)明顯的例子。而且,由于這種欠采樣,曾困擾Hubble望遠(yuǎn)鏡高精度天體測量工作長達(dá)10年之久。Anderson和King提出了有效點(diǎn)擴(kuò)散函數(shù)(ePSF-effective Point——Spread Function)新概念[6],并利用抖動(dòng)的圖像序列提取精確的ePSF,基本上解決了Hubble望遠(yuǎn)鏡圖像欠采樣的問題。然而,Anderson本人也承認(rèn)求解ePSF的復(fù)雜性。當(dāng)圖像足夠采樣時(shí),ePSF擬合法與經(jīng)典的高斯擬合法相比,不會(huì)有精度上的增益。
為了更深刻地認(rèn)識(shí)像素相位誤差對(duì)天體測量的影響,在生成仿真圖像上利用高斯函數(shù)作為點(diǎn)擴(kuò)散函數(shù)模型,探究了影響像素相位誤差的決定因素,也分析了像素相位誤差對(duì)位置測量精度的影響。
像素相位是星像中心相對(duì)于像素邊界的位置,方向的像素相位可分別定義為φ(x)=x-int(x+0.5),φ(y)=y-int(y+0.5),(x、y分別為星像中心的橫坐標(biāo),縱坐標(biāo))。采用文 [6]中的圖來說明像素相位誤差是如何產(chǎn)生的。圖1顯示的是兩個(gè)不同的點(diǎn)擴(kuò)散函數(shù)(PSF)模型擬合相同的一維欠采樣星像輪廓,這個(gè)星像輪廓由一個(gè)直方圖表示,其中包含星像中心的3 pix。圖中可以看出同樣精確擬合星像輪廓的兩個(gè)不同的點(diǎn)擴(kuò)散函數(shù)模型擬合出的星像中心有明顯差異,達(dá)到0.07 pix,這個(gè)偏差就是像素相位誤差。

圖1 這個(gè)直方圖顯示欠采樣的一維星像輪廓最中心3 pix的灰度。實(shí)線表示純粹的高斯模型擬合結(jié)果,虛線表示較為尖銳的高斯模型擬合結(jié)果,并且擬合的范圍較大。箭頭顯示了這兩個(gè)高斯模型在峰值的位置上相差0.07 pix(原圖來自文[6])Fig.1 Flux distribution int the x-direction cross section through the center of a stellar profile.The histogram shows the pixel values for the innermost 3 pixels of the underampled one-dimensional stellar profile.The solid curve is the fit with a pure Gaussian model.The dotted curve is the fit with a sharper Gaussian for a broader range(adopted from Ref[6])
為解決Hubble空間圖像欠采樣的問題,美國學(xué)者Anderson和King提出了ePSF的概念[6]。假如一個(gè)星像中心位于(xc,yc),(i,j)是該星像周圍的一個(gè)像素,它們與ePSF存在如下關(guān)系,

其中pij為像素(i,j)的灰度值;S*為天空背景值;f*為星像的通量因子; ψE(Δx,Δy)=ψE(i-xc,jyc)為點(diǎn)(i,j)相對(duì)于星像中心的有效點(diǎn)擴(kuò)散函數(shù)值(ePSF)。
(1)式的提出給仿真圖像的生成帶來了很大方便,不必像Stone[1]那樣在一個(gè)像素范圍內(nèi)進(jìn)行積分求得該像素的灰度值,而只需要利用(1)式進(jìn)行簡單的計(jì)算就可以生成像素點(diǎn) (i,j)處的灰度值。本文正是基于ePSF的概念,生成所有的仿真圖像。
用高斯函數(shù)(Gaussian)模型描述星像輪廓:

I(x,y)為像素(x,y)處的灰度值;H是高斯函數(shù)的峰值;B為天空背景;(xc,yc)是星像中心的坐標(biāo)位置;R為高斯函數(shù)的標(biāo)準(zhǔn)差,它與大氣視寧度有關(guān)。R與星像的半高全寬(FWHM)存在如下關(guān)系:

比較(1)式和(3)式發(fā)現(xiàn),ePSF和Gaussian存在如下關(guān)系:

對(duì)仿真圖像,采用2k×2k的分辨率。生成仿真圖像的主要步驟如下,首先設(shè)置圖像的高度、寬度及調(diào)用CFITSIO庫[7]創(chuàng)建背景圖像。然后在背景圖像中通過設(shè)置星像的中心位置、通量、半高全寬來添加星像,其中星像采用高斯模型(式2)。最后對(duì)圖像添加泊松噪聲。圖2是仿真圖像的一部分。
這里生成的星像半高全寬都為1.5 pix的50幅仿真圖像,其中每幅圖像添加像素相位隨機(jī)的200顆亮星。采用亮星的目的是亮星具有很高的信噪比,可以更好地觀察像素相位誤差的影響。
現(xiàn)在采用高斯函數(shù)(式2)來實(shí)驗(yàn)分析什么是影響像素相位誤差的決定因素及怎樣獲得精確的高斯函數(shù)模型。高斯函數(shù)擬合過程中,如果標(biāo)準(zhǔn)差R不參與迭代,R為固定值,大小為設(shè)置的初值,這時(shí)擬合只需求解B、H、xc、yc4個(gè)參數(shù)。如果R參與迭代,則求解的參數(shù)有5個(gè),分別是B、H、xc、yc、R。

圖2 生成仿真圖像的一部分Fig.2 A part of the simulated image
迭代擬合求4個(gè)參數(shù)B、H、xc、yc的解。它們的初始值分別設(shè)為B0、H0、x0、y0,其中x0、y0是目標(biāo)框(圈定星像的矩形框)內(nèi)最大灰度值的坐標(biāo)位置,B0是目標(biāo)框內(nèi)的最小灰度值,H0是目標(biāo)框內(nèi)最大灰度值與最小灰度值之差。在迭代過程中有下述關(guān)系:xc=xc+Δx,yc=yc+Δy,B=B+ΔB,H=H+ΔH。等式中所求得的參數(shù)值只是某次迭代過程中的值,并不是各個(gè)參數(shù)的確切值。要得到較為準(zhǔn)確的值,需要多次迭代。迭代求解結(jié)束的條件是: |Δx|<0.001,|Δy|<0.001,|ΔB|<0.01,|ΔH|<0.01。在實(shí)驗(yàn)過程中對(duì)高斯函數(shù)標(biāo)準(zhǔn)差R分別設(shè)置5個(gè)初值,R1=FWHM0/2.40(FWHM0為星像的半高全寬1.5)、R2=FWHM0/2.38、R3=FWHM0/2.35、R4=FWHM0/2.33、R5=FWHM0/2.30。
針對(duì)上述R的不同取值,通過觀察這50幅仿真圖像計(jì)算出每個(gè)星像的殘差(xc-x0,yc-y0),每幅圖像200顆星,即總共有10000個(gè)殘差,用這些殘差描述像素相位誤差。用圖3~7表示像素相位誤差相對(duì)于像素相位的變化,其中橫軸是x、y方向的像素相位,縱軸是對(duì)應(yīng)方向的像素相位誤差。

圖3 R1=FWHM0/2.40,模型的半高全寬相對(duì)于星像真實(shí)的半高全寬偏小,殘差總體上先負(fù)后正Fig.3 Fitting residuals of stellar positions for R1=FWHM0/2.40,under which the model FWHM are less than the true FWHM.The residuals overall have the same signs as pixel phases

圖4 R2=FWHM0/2.38,模型的半高全寬相對(duì)于星像真實(shí)的半高全寬也是偏小Fig.4 Fitting residuals of stellar positions for R2=FWHM0/2.38,under which the model FWHM are again less than the true FWHM

圖5 R3=FWHM0/2.35模型的半高全寬與星像真實(shí)的半高全寬最接近,殘差最小Fig.5 Fitting residuals of stellar positions for R3=FWHM0/2.35,under which the model FWHM are closest to the true FWHM and the residuals reach minima

圖6 R4=FWHM0/2.33模型的半高全寬相對(duì)于星像真實(shí)的半高全寬偏大,殘差總體上先正后負(fù)Fig.6 Fitting residuals of stellar positions for R4=FWHM0/2.33,under which the model FWHM are larger than the true FWHM.The residuals overall have an anti-correlation with the pixel phases

圖7 R5=FWHM0/2.30模型的半高全寬相對(duì)于星像真實(shí)的半高全寬也是偏大Fig.7 Fitting residuals of stellar positions for R5=FWHM0/2.30,under which the model FWHM are again larger than the true FWHM
實(shí)驗(yàn)發(fā)現(xiàn),模型的半高全寬是影響像素相位誤差的決定因素,當(dāng)模型的半高全寬越接近星像真實(shí)的半高全寬時(shí),像素相位誤差越小。圖5顯示出高斯模型的半高全寬與星像真實(shí)的半高全寬偏差最小時(shí),像素相位誤差也最小。把圖3~7的像素相位誤差幅度的平均曲線見圖8。其中R1、R5分別表示高斯模型的半高全寬相對(duì)于星像真實(shí)的半高全寬偏小與偏大的情況,它們對(duì)應(yīng)曲線的幅度很大且相位相差180°。
另外還發(fā)現(xiàn)x方向的像素相位誤差只受x方向的像素相位影響,而與y方向的像素相位沒有關(guān)系,反之亦然。取R5=FWHM0/2.30的情況進(jìn)行實(shí)驗(yàn)發(fā)現(xiàn),y方向的像素相位誤差對(duì)于x方向的像素相位是隨機(jī)分布的(圖9左),x方向的像素相位誤差對(duì)于y方向的像素相位也是隨機(jī)分布的(圖9右)。

圖8 橫軸為像素相位,縱軸為像素相位誤差的平均值。其中R1、R5的情況是模型的半高全寬相對(duì)星像真實(shí)的半高全寬偏得較大,此時(shí)像素相位誤差的幅度也很大。R3是幾個(gè)值中偏得最小的,此時(shí)像素相位誤差幅度也最小Fig.8 Average pixel-phase error versus pixel phase in the x direction.The R1and R5cases are those large with large deviations of the model FWHM from the true FWHM and amplitudes of pixel-phase error variations.In the R3case,the model FWHM are closest to the true FWHM and the amplitude of pixel-phase error variation reaches minimum
此時(shí)迭代擬合要求解B、H、xc、yc、R這5個(gè)參數(shù)。與R不參與迭代情況唯一不同的是R的最終解是進(jìn)行多次迭代后的結(jié)果,而不僅僅是設(shè)置初值。設(shè)R的初始值為FWHM0/2.30,且R退出迭代的條件是<0.001。每迭代一次,高斯函數(shù)的標(biāo)準(zhǔn)差就趨近于星像真實(shí)的標(biāo)準(zhǔn)差FWHM0/,經(jīng)過幾次迭代之后,兩者大小幾乎相等,也就是迭代后的高斯模型半高全寬幾乎和星像真實(shí)的半高全寬相等。實(shí)驗(yàn)發(fā)現(xiàn),R參與迭代的高斯擬合效果非常好,殘差很小(圖10)。可見,擬合模型精確與否取決于模型的半高全寬與星像真實(shí)的半高全寬是否趨于相等。

圖9 (a)y方向的殘差關(guān)于x方向像素相位的函數(shù)(b)x方向的殘差關(guān)于y方向像素相位的函數(shù)Fig.9 (a)Residuals of y positions versus x pixel phases(b)Residuals of x positions versus y pixel phases

圖10 參與迭代的高斯擬合效果(a)x方向的殘差關(guān)于x方向像素相位的函數(shù) (b)y方向的殘差關(guān)于y方向像素相位的函數(shù)Fig.10 Results of Gaussian fitting with our iterative approach involving the R parameter.(a)Residuals of x positions versus x pixel phases (b)Residuals of y positions versus y pixel phases
和上節(jié)生成圖像的情況不一樣,生成星像的半高全寬分別是 1.4、1.5、1.6、1.7、1.9、2.0、2.3、2.5、2.7、3.0的一系列仿真圖像。采用不精確的高斯函數(shù)模型(標(biāo)準(zhǔn)差R不參與迭代,其大小為FWHM0/2.30)來實(shí)驗(yàn)分析不同半高全寬下像素相位誤差幅度的大小。實(shí)驗(yàn)發(fā)現(xiàn),像素相位誤差的幅度隨著星像半高全寬的增大而變小,當(dāng)半高全寬大于2 pix,像素相位誤差的幅度幾乎為0(圖11)??梢?,在星像的半高全寬大于2 pix,高斯擬合法無論標(biāo)準(zhǔn)差精確與否都不再受像素相位誤差的影響。

圖11 橫軸是星像的半高全寬,縱軸是像素相位誤差的幅度Fig.11 Amplitude of pixel-phase error variation versus true FWHM
生成的星像的半高全寬為1.4~2.2,且兩相鄰值間距為0.1的一系列仿真圖像,其中每幅圖像添加位置和通量都隨機(jī)的星(更符合實(shí)際圖像的特征)。同樣采用高斯擬合法,并在實(shí)驗(yàn)過程中分3種情況(R參與迭代,R不參與迭代且大小為FWHM /,R不參與迭代且大小為FWHM/2.30)進(jìn)行實(shí)驗(yàn),比較這3種方法在不同半高全寬下的位置測量精度。先采用李展[8]等人的方法歸算出每一顆星像位置的標(biāo)準(zhǔn)差,為了方便,進(jìn)一步用所有星像位置的平均標(biāo)準(zhǔn)差來反映算法的精度。圖12表示這3種方法擬合星像的位置精度。其中橫軸為星像的半高全寬,左圖和右圖的縱軸分別代表x,y方向擬合星像位置的平均標(biāo)準(zhǔn)差大小,顯然平均標(biāo)準(zhǔn)差越小,精度越高。
實(shí)驗(yàn)發(fā)現(xiàn):(1)R不參與迭代其大小取FWHM/2.30時(shí),高斯模型的半高全寬和星像真實(shí)的半高全寬偏得較大,此時(shí)高斯模型是不精確的。在半高全寬小于2 pix時(shí)會(huì)受像素相位誤差的影響而導(dǎo)致此方法的位置測量精度不高。(2)R不參與迭代其大小取FWHM /時(shí),高斯模型的半高全寬等于星像真實(shí)的半高全寬,此時(shí)認(rèn)為高斯模型等價(jià)于ePSF模型(式4),這種情況沒有像素相位誤差。此時(shí)R參與迭代的高斯模型與ePSF模型都是精確的擬合模型,它們在各半高全寬下精度近似相等。(3)隨著半高全寬的增大,像素相位誤差的影響程度越來越小,在星像的半高全寬大于2 pix時(shí),3種方法的精度近似相等(圖12)。

圖12 (a)x方向的平均標(biāo)準(zhǔn)差大小 (b)y方向的平均標(biāo)準(zhǔn)差大小Fig.12 (a)Standard deviations of x positions versus true FWHM (b)Standard deviations of y positions versus true FWHM
為了深入探究像素相位誤差的規(guī)律,生成一系列仿真圖像并采用二維高斯函數(shù)模型作為點(diǎn)擴(kuò)散函數(shù)模型。實(shí)驗(yàn)表明:(1)擬合模型的半高全寬與星像真實(shí)的半高全寬偏差越大,像素相位誤差的幅度越大,反之亦真。(2)當(dāng)星像真實(shí)的半高全寬大于2 pix時(shí),像素相位誤差可以忽略不計(jì)。(3)當(dāng)高斯模型的半高全寬參與迭代擬合時(shí),高斯與ePSF這兩個(gè)函數(shù)模型擬合的位置精度近似相等。
致謝:感謝暨南大學(xué)計(jì)算機(jī)系孟小華、李展老師在實(shí)驗(yàn)中提出的建設(shè)性意見和提供的幫助。
[1]Stone R C.A Comparison of Digital Centering Algorithms[J].AJ,1989(97):1227 - 1237.
[2]Zacharias N,Urban S E,Zacharias M I,et al.The First Us Naval Observatory CCD Astrograph Catalog(UCAC1) [J].AJ,2000,120(4):2131 - 2147.
[3]Zacharias N,Urban S E,Zacharias M I,et al.The Second Us Naval Observatory CCD Astrograph Catalog(UCAC2) [J].AJ,2004,127(5):3043 - 3059.
[4]Zacharias N,F(xiàn)inch C,Girard T,et al.The Third US Naval Observatory CCD Astrograph Catalog(UCAC3) [J].AJ,2010,139(6):2184 - 2199.
[5]Zacharias N.UCAC3 Pixel Processing [J].AJ,2010,139(6):2208 - 2217.
[6]Anderson J,King I R.Toward High Precision Astrometry with WFPC2.I.Deriving an Accurate Point-Spread Function [J].The Publications of the Astronomical Society of the Pacific,2000,112(776):1360-1392.
[7]CFITSIO [EB/OL].http://heasarc.gsfc.nasa.gov/.
[8]李展,彭青玉,韓國強(qiáng).CCD圖像數(shù)字定心算法的比較 [J].天文學(xué)報(bào),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.