張 鑫

(1.中國科學(xué)院國家天文臺空間天文與技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100101；2.中國科學(xué)院大學(xué)，北京 100049)

天文圖像模擬是將圖像處理技術(shù)與天文學(xué)模型相結(jié)合，這種模擬在天文學(xué)研究中已經(jīng)越來越廣泛使用。模擬有遇見未來之功效，通過模擬可以對未來的設(shè)備進(jìn)行更加科學(xué)的評估?；谶@種原因，越來越多的人開始從事天文圖像模擬的研究。在對天文的觀測研究中，有一些被學(xué)者和觀測者熟識的天文觀測模擬軟件，如WorldWide Telescope，Astronomy CCD Calculator，Stellarium等，這些軟件對天體空間位置、天體特性的描述非常細(xì)致，給觀測者帶來了很大的便利，但是對于遙遠(yuǎn)星系的描述這些軟件還不是很完善。隨著大型觀測設(shè)備的出現(xiàn)，人們對未知天空的觀測變得越來越遙遠(yuǎn)，星系模擬也成為評估大型天文設(shè)備很重要的一部分。文[1]介紹了SkyMaker，該軟件很細(xì)致地描述了點(diǎn)擴(kuò)散函數(shù)，包括大氣擾動、望遠(yuǎn)鏡運(yùn)動、像差、衍射等各種因素，星系形態(tài)特征利用了De Vaucouleurs方法進(jìn)行描述。文[2]利用Shapelet的方法對星系的形態(tài)進(jìn)行模擬，該方法對星系形態(tài)采用多維模擬，文中與小波(Wavelet)等方法進(jìn)行比較，對殘差進(jìn)行分析，證明其方法的可靠性。LSST[3]成立了專門的圖像模擬團(tuán)隊(duì)，該團(tuán)隊(duì)主要工作是為LSST做模擬。他們的模擬針對全天區(qū)，利用宇宙模型建立星表確定天體的分布，并且建立了多種天體現(xiàn)象的模型(天體模型包括恒星、星系、小行星等)，在傳輸過程中考慮天氣情況、大氣擾動等，最后還考慮了結(jié)合巡天設(shè)備影響產(chǎn)生的點(diǎn)擴(kuò)散函數(shù)。利用準(zhǔn)確的模擬結(jié)果可以對LSST系統(tǒng)的各項(xiàng)性能進(jìn)行分析處理。另外，LSST對天體形態(tài)的模擬采用蒙特卡羅的方法，模擬光子的射入。文[4]提到了模擬的兩種不同方法，一種是類似LSST模擬光子射入的方法，另一種是按照像素模擬的方法，并比較兩種方法的差異，文中用了簡單的數(shù)學(xué)模型，主要強(qiáng)調(diào)處理的速度。

本文的模擬利用哈勃(Hubble)望遠(yuǎn)鏡的數(shù)據(jù)和觀測條件。圖像模擬工作包含：(1)模擬天體的分布和天體的形態(tài)，主要是星系的形態(tài)；(2)模擬觀測條件以及各種噪聲，包括天光背景、暗流、讀出噪聲等；(3)模擬儀器的點(diǎn)擴(kuò)散函數(shù)；(4)對觀測結(jié)果進(jìn)行統(tǒng)計(jì)分析。圖1為圖像模擬的流程圖。

本文的主要目標(biāo)是利用簡單通用的模型對天文圖像實(shí)現(xiàn)快速的模擬。在程序編寫上應(yīng)用更多的優(yōu)化，使程序達(dá)到快速可用的目的。本文還利用SExtractor對模擬結(jié)果進(jìn)行提取，將提取的結(jié)果與參考結(jié)果進(jìn)行比較分析，通過分析確定模擬結(jié)果在實(shí)際應(yīng)用上的可行性。

1 圖像模擬的關(guān)鍵模型及方法

本文模擬利用HUDF ACS WFC i(F775)波段的數(shù)據(jù)，該數(shù)據(jù)生成的圖像約為10 000×10 000像素，由于圖像在空間上存在一定的旋轉(zhuǎn)，圖像4個角包含未拍攝的區(qū)域，為了實(shí)驗(yàn)方便，截取了其中4 500×4 500像素的中間區(qū)域，面積約為5.1 arcmin2。采用天文學(xué)中一些簡單高效的模型，主要為了對數(shù)據(jù)進(jìn)行快速處理。在觀測條件上，模擬了哈勃望遠(yuǎn)鏡在600 km高空的觀測條件。

1.1 星系的模型

星系的模型采用Sersic模型[5]，數(shù)學(xué)表達(dá)式如下：

其中，I(R)為R處的流量密度；Ie為有效半徑處的流量密度；Re為有效半徑；n為Sersic系數(shù)，通過該系數(shù)可以控制星系的彌散程度；bn為常數(shù)。通過變換可以將(1)式變換成如下的表達(dá)式：

其中，I0為中心位置的流量密度；n為Sersic系數(shù)；a為常數(shù)，表示當(dāng)R為a時的流量密度為中心位置流量密度的e-1倍。

圖1 圖像模擬流程圖Fig.1 Flow Diagram of Image Simulation

若星系為圓形區(qū)域，對(2)式在圓內(nèi)積分，得到星系的流量：

其中，Γ為gamma函數(shù)。

在模擬中，通常用橢率表示星系的形狀，這就需要將星系的流量密度分布在整個橢圓平面內(nèi)，定義橢圓的橢率為e，則(2)式中R可以用歸一化的橢圓半徑表示：

其中，x，y表示坐標(biāo)，坐標(biāo)的原點(diǎn)為橢圓的中心。

星系在圖像中的表現(xiàn)，除了彌散程度和形狀的特性外，還存在方向性，這就需要對已經(jīng)計(jì)算出來的橢圓形的Sersic分布進(jìn)行坐標(biāo)變換，假設(shè)星系坐標(biāo)系(長軸為x軸，短軸為y軸)與圖像坐標(biāo)系的夾角為θ，則變換矩陣為

新的坐標(biāo)位置(x′， y′)為

通過(1)～(6)式的計(jì)算可以模擬星系的形狀、大小、方向等要素，而對星系的模擬需要的參數(shù)，如n(Sersic系數(shù))、Re(有效半徑)、亮度等都是通過對哈勃超深空場ACS WFC的數(shù)據(jù)進(jìn)行提取獲得。在文[6]中利用SExtractor對哈勃超深空場的數(shù)據(jù)進(jìn)行提取，并和文[7]提供的數(shù)據(jù)進(jìn)行對比，結(jié)果比較理想，文中用到的參數(shù)均來自于文[6]提取的星表。

1.2 觀測條件及噪聲的模擬

對于觀測條件的模擬，在圖像中的表現(xiàn)其實(shí)就是對噪聲的還原，還原的真實(shí)與否是影響模擬圖像效果的最主要因素。觀測產(chǎn)生的噪聲主要包括固有噪聲、信號噪聲等。

固有噪聲是指源于探測器自身的噪聲，包括熱噪聲、散粒噪聲、生成-復(fù)合噪聲和閃變噪聲，在文中用讀出噪聲表示。為了和哈勃的數(shù)據(jù)進(jìn)行對比，讀出噪聲利用哈勃的曝光時間計(jì)算器[8]計(jì)算得出。

信號噪聲主要包含背景噪聲和由于信號的量子性產(chǎn)生的隨機(jī)噪聲。背景噪聲利用了哈勃望遠(yuǎn)鏡測得的天光背景噪聲，包括黃道光和赤道光。

對于信號，由于其量子性，到達(dá)探測器會有一定的隨機(jī)性，基于這種特性也會產(chǎn)生一定的噪聲，在文中用泊松分布表示。

1.3 點(diǎn)擴(kuò)散函數(shù)模擬

對于地面望遠(yuǎn)鏡，影響點(diǎn)擴(kuò)散函數(shù)的主要因素有大氣擾動和儀器自身的光學(xué)特性。而對于哈勃望遠(yuǎn)鏡這樣的空間設(shè)備，觀測條件要遠(yuǎn)遠(yuǎn)優(yōu)于地面望遠(yuǎn)鏡，因此對點(diǎn)擴(kuò)散函數(shù)的影響僅僅來自觀測設(shè)備的光學(xué)特性。

文中的模擬是對空間望遠(yuǎn)鏡成像模擬，點(diǎn)擴(kuò)散函數(shù)采用高斯(Gaussian)函數(shù)表示：

其中，I0為常數(shù)，表示中心位置的流量密度，可以控制高斯函數(shù)的幅度；[x，y]為圖像坐標(biāo)；σ用來控制圖像的質(zhì)量。對于高斯函數(shù)通?？梢酝ㄟ^所給的半高全寬(FWHM)求得：

對于地面望遠(yuǎn)鏡通常需要考慮大氣擾動的效應(yīng)，所以模擬中用到的點(diǎn)擴(kuò)散函數(shù)需要考慮觀測設(shè)備和大氣擾動的綜合效應(yīng)，對于這類點(diǎn)擴(kuò)散函數(shù)模擬通常采用Moffat[9]模型：

其中，I0為中心位置處的流量密度；α和β為尺度因子，可以控制函數(shù)的彌散程度；[x，y]為圖像坐標(biāo)。α和β關(guān)系如下：

文[10]中提到，當(dāng)β＝4.765時，Moffat模型可以很好地模擬大氣擾動。

1.4 子像素級處理

模擬的目的之一是對未來某一個觀測設(shè)備的評估，文中模擬利用的數(shù)據(jù)只符合Hubble ACS的特性，例如像素的大小為0.03×0.03 arcsec2，為了使程序能夠?yàn)槠渌^測設(shè)備評估利用，模擬需要做到子像素級。子像素級處理主要應(yīng)用在天體形態(tài)、亮度的模擬過程中，根據(jù)需要將原像素分成若干個子像素，根據(jù)天體的特性，在子像素中計(jì)算得到該像素的值，計(jì)算量十分龐大，這就需要對程序進(jìn)行優(yōu)化。

1.5 并行計(jì)算

在計(jì)算天體的形態(tài)、亮度等要素時，由于將原有像素拆分為子像素計(jì)算，這樣導(dǎo)致計(jì)算量倍增。另外，每一個天體都是獨(dú)立的個體，單獨(dú)計(jì)算對其他天體不會產(chǎn)生任何影響。所以，文中將計(jì)算每一個天體的形態(tài)、亮度等要素進(jìn)行獨(dú)立計(jì)算，這樣，將模擬中的所有天體按照中央處理器的數(shù)目進(jìn)行平均分配，不同的中央處理器計(jì)算的結(jié)果返給主進(jìn)程進(jìn)行合并，形成最終模擬的整幅圖像。

2 結(jié)果分析

2.1 參考樣本分析

利用文[6]提供的數(shù)據(jù)，其中包含了提取的Sersic系數(shù)，而Beckwith提供了一份更加準(zhǔn)確的數(shù)據(jù)，不過在該數(shù)據(jù)中對于天體彌散程度的描述只給出了高斯參數(shù)，但是其位置和亮度更加準(zhǔn)確。事實(shí)上，文[6]的數(shù)據(jù)也給出了一些與Bechwith數(shù)據(jù)的對比，例如，前者的數(shù)據(jù)給出了中心位置與Bechwith的對比，本文對兩個數(shù)據(jù)做了更詳細(xì)的對比。

圖2是對兩份數(shù)據(jù)提取天體亮度的對比，藍(lán)色線為文[6]數(shù)據(jù)天體亮度的統(tǒng)計(jì)，紅色線為Bechwith的數(shù)據(jù)天體亮度的統(tǒng)計(jì)。從圖2可以看出，兩份數(shù)據(jù)的亮度存在一定偏差，文[6]的數(shù)據(jù)天體亮度要比實(shí)際偏大。

圖3是對兩個樣本有效半徑的對比結(jié)果，橫坐標(biāo)是提取天體的亮度信息，縱坐標(biāo)是Benchwich的數(shù)據(jù)中天體的有效半徑與文[6]的數(shù)據(jù)中相應(yīng)的天體有效半徑之差。從圖3可以看出，差值絕大部分集中在0值附近，但是文[6]的數(shù)據(jù)有部分半徑偏大。

圖2 Coe D.等提供的數(shù)據(jù)與Bechwith的數(shù)據(jù)天體亮度對比Fig.2 AB magnitude comparison between Coe et al.catalog and Bechwith et al.catalog

圖3 兩種參考樣本提取天體有效半徑的對比Fig.3 Effective Radius comparison between Coe et al.catalog and Bechwith et al.catalog

綜合上述對樣本的對比，利用文[6]的數(shù)據(jù)模擬的圖像與真實(shí)圖像存在亮度偏差，另外，對于天體的大小也有部分天體與真實(shí)圖像存在一定差別。從這一點(diǎn)，也可以解釋后面章節(jié)中實(shí)驗(yàn)結(jié)果圖4的模擬結(jié)果存在一些比原圖大的天體的現(xiàn)象。

2.2 模擬結(jié)果及分析

圖4展示了模擬的結(jié)果，圖4(a)是模擬圖像，圖4(b)是經(jīng)過MultiDrizzle程序處理后的ACS F775波段的圖像。由于原圖在空間存在一定角度，并且4個角上沒有目標(biāo)信息，所以本文用到的數(shù)據(jù)和圖像是對原圖的截取，截取了原圖x方向3001到7500像素，y方向3001到7500像素的方形區(qū)域，總共4 500×4 500個像素，經(jīng)過MultiDrizzle程序處理后的ACS F775波段的圖像每個像素的視場約為0.03×0.03 arcsec2，模擬區(qū)域的視場為5.1 arcmin2。哈勃望遠(yuǎn)鏡團(tuán)隊(duì)給出ACS WFC i(F775)波段零點(diǎn)星等(zeropoint)的值為25.654，這里零點(diǎn)的意思是流量為1e-時天體的亮度值。根據(jù)(11)式可以計(jì)算天體的流量：

由(11)式計(jì)算的流量值是信號的流量。而噪聲包括天光背景、暗流、讀出噪聲等，是利用哈勃望遠(yuǎn)鏡的曝光時間計(jì)算器計(jì)算得到。計(jì)算結(jié)果以e-為單位的流量。所有計(jì)算的結(jié)果除以增益(單位，e-/DN)得到表現(xiàn)在圖像上的數(shù)字值。

圖4 模擬圖像結(jié)果Fig.4 Simulation Image vs.Original Image(a： Simulation Image， b： Original Image)

HUDF ACS WFC I波段的極限星等為29等，所以本文提取的截止亮度到29等。文中用SExtractor作為提取工具，SExtractor的提取主要參數(shù)見附錄。表1是對提取數(shù)目的統(tǒng)計(jì)，參考的對象為文[6]的數(shù)據(jù)。表中有4項(xiàng)數(shù)據(jù)，第1項(xiàng)是對原圖提取數(shù)據(jù)的統(tǒng)計(jì)，第2項(xiàng)是對模擬圖像提取數(shù)目的統(tǒng)計(jì)，第3項(xiàng)是粗略地做了匹配統(tǒng)計(jì)，匹配的方法是將模擬數(shù)據(jù)提取的天體中心與文[6]數(shù)據(jù)的天體中心作比較，選出距離最接近的一對被認(rèn)為是同一個天體。第4項(xiàng)是計(jì)算匹配率的比例，該比例通過(12)式計(jì)算得到：

其中，rmatch為匹配率，Nref為參考數(shù)據(jù)提取天體的數(shù)目；Nmatch為對模擬數(shù)據(jù)提取天體的數(shù)目。從表1匹配的數(shù)目可以看出，這種針對位置的匹配度均在85%以上，考慮到噪聲影響和提取方法的差異，這種匹配程度是可以接受的。

表1 利用模擬圖像提取的天體與HUDF catalog天體匹配數(shù)目統(tǒng)計(jì)Table 1 Extracted galaxies analysis

圖5是對小于29等的天體數(shù)目的統(tǒng)計(jì)，從圖5可以看出，對于模擬圖像提取的天體數(shù)目和參考數(shù)據(jù)天體的數(shù)目隨亮度變化的趨勢是比較一致的。

圖6是將模擬數(shù)據(jù)提取的天體有效半徑與參考數(shù)據(jù)天體的有效半徑進(jìn)行對比，采用了計(jì)算匹配天體有效半徑之差的方法進(jìn)行對比，匹配天體的方法見表1，在這里只是確定了位置上的匹配，這兩個統(tǒng)計(jì)的樣本都是亮度小于29等的天體。圖6(a)統(tǒng)計(jì)該差值隨著亮度分布的變化，橫坐標(biāo)是亮度，縱坐標(biāo)是參考數(shù)據(jù)中天體的有效半徑與相應(yīng)的模擬數(shù)據(jù)提取的天體有效半徑之差，從分布圖上可以看出，亮度越小，有效半徑的差值越小。圖6(b)是對差值的分布進(jìn)行統(tǒng)計(jì)，從統(tǒng)計(jì)可以看出，這一差值絕大多數(shù)分布在0值左右，有部分存在較大的偏差。導(dǎo)致這種差異的原因有很多，其中包括提取樣本的精度(2.1節(jié)中對兩種參考樣本的分析)，圖像模擬中加入的信號噪聲、背景噪聲和設(shè)備噪聲及利用SExtractor對天體進(jìn)行提取的參數(shù)也存在一定差異。

圖7是對模擬圖像提取的天體亮度與參考數(shù)據(jù)的天體亮度進(jìn)行對比，對比的方法同圖6(a)，通過亮度之差隨天體亮度的分布進(jìn)行分析。圖中，橫坐標(biāo)是亮度，縱坐標(biāo)是參考數(shù)據(jù)中天體的亮度與相應(yīng)的模擬數(shù)據(jù)中提取的天體亮度之差。統(tǒng)計(jì)的樣本依然是亮度小于29等天體。從圖7可以看出，這一亮度差值大部分分布在0軸上下，偏差絕大多數(shù)集中在0.5以下。從圖7還可以看出一個規(guī)律，對于亮度越大的天體，模擬后再提取，其亮度值越準(zhǔn)確。

圖5 提取天體數(shù)目統(tǒng)計(jì)對比Fig.5 Extracted Galaxies from Simulation Image vs.HUDF catalog

圖6 模擬圖像所提取天體的有效半徑分析圖Fig.6 Analysis of Effective Radius for Galaxies from Simulation Image

圖7 模擬圖像中提取天體亮度分析圖Fig.7 Magnitude analysis for Galaxies from Simulation Image

2.3 性能分析

本文中圖像模擬軟件以及統(tǒng)計(jì)的數(shù)據(jù)生成都是用JAVA程序開發(fā)。在天文中用JAVA開發(fā)的程序較少，本文之所以選擇JAVA程序開發(fā)主要因?yàn)镴AVA程序存在以下優(yōu)勢：(1)JAVA去掉了指針，有自動的垃圾回收機(jī)制，這樣減小了編寫程序時由于內(nèi)存泄露出現(xiàn)的問題，因此，能夠?qū)⒏嗟臅r間用在算法的編寫上；(2)JAVA程序是面向?qū)ο蟮木幊?，在程序的可讀性上要好于C程序。(3)JAVA支持并行計(jì)算，且編寫簡單，對于本程序的并行計(jì)算有很大的幫助。

另外，本文的程序運(yùn)行在Linux系統(tǒng)下，硬件配置CPU為Intel core i5 4核3.1 GHz，內(nèi)存為16 GB，在程序運(yùn)行時，4核全部使用，計(jì)算天體的數(shù)目為4 107，計(jì)算所需要的時間為1 267 s，平均每個天體所需要時間為308 ms。程序運(yùn)行時峰值所需要的內(nèi)存為4.1 GB。

程序運(yùn)行的大部分時間和內(nèi)存都是消耗在利用Sersic模型計(jì)算天體的形態(tài)，在計(jì)算天體的形態(tài)時，由于1.4節(jié)所述子像素的需求以及計(jì)算的精度，需要將一個像素分解成若干個小像素，計(jì)算時間和內(nèi)存的消耗與像素分解的個數(shù)成正比。

對于星系形態(tài)的計(jì)算采用了二維的數(shù)值積分，為了計(jì)算更加準(zhǔn)確，本文將一個像元分解為若干個子像元，通過利用子像元積分得到每一個像元的準(zhǔn)確數(shù)據(jù)，文中兼顧性能需求沒有對像元進(jìn)行無限劃分，只是分解為32×32個子像素，4 500×4 500像素的圖像就變成了144 000×144 000像素的圖像，這樣數(shù)據(jù)量已經(jīng)很大，采用并行及算法上的優(yōu)化使得每個天體的計(jì)算時間達(dá)到308 ms。

3 總結(jié)及展望

本文對圖像模擬用到的方法都是天文中比較基礎(chǔ)的模型，例如Sersic模型、高斯模型。將這些基礎(chǔ)模型應(yīng)用到模擬中，可以達(dá)到簡單高效的目的。但是簡單高效只是在模擬結(jié)果真實(shí)正確的基礎(chǔ)上才能成立。通過本文第2部分的分析可以發(fā)現(xiàn)，模擬的結(jié)果與真實(shí)數(shù)據(jù)存在一定差異，但是基本特征都能夠正確地表現(xiàn)(例如大小、方向、亮度等)；通過統(tǒng)計(jì)分析可以看出，對于模擬圖像提取的天體的數(shù)目、亮度、大小與參考數(shù)據(jù)比較，錯誤的數(shù)量都在可以控制的范圍內(nèi)。另外，本文對模擬程序?qū)崿F(xiàn)了并行計(jì)算，而且在程序中已做了一些優(yōu)化(例如，三角函數(shù)用插值方法計(jì)算)，通過對性能的分析可以看到，計(jì)算天體的平均時間達(dá)到了毫秒級，再次說明了本文的模擬成像速度快。

模擬程序中，點(diǎn)擴(kuò)散函數(shù)只用到高斯函數(shù)，雖然簡單，但是對于反映圖像的真實(shí)性還存在一定的差別。另外，模擬的目的主要是為一些新的巡天項(xiàng)目分析服務(wù)，所以下階段的目標(biāo)是在模擬程序上對點(diǎn)擴(kuò)散函數(shù)進(jìn)行改進(jìn)，加入真實(shí)的光學(xué)系統(tǒng)特性。另外需要引入光譜的模擬以適應(yīng)多波段成像模擬。同時，在分析上，將模擬與新的巡天項(xiàng)目相結(jié)合，得到與不同硬件相關(guān)的分析。