一种快速的天文图像模拟方法及分析*

张 鑫

(1.中国科学院国家天文台空间天文与技术重点实验室，北京 100101；2.中国科学院大学，北京 100049)

天文图像模拟是将图像处理技术与天文学模型相结合，这种模拟在天文学研究中已经越来越广泛使用。模拟有遇见未来之功效，通过模拟可以对未来的设备进行更加科学的评估。基于这种原因，越来越多的人开始从事天文图像模拟的研究。在对天文的观测研究中，有一些被学者和观测者熟识的天文观测模拟软件，如WorldWide Telescope，Astronomy CCD Calculator，Stellarium等，这些软件对天体空间位置、天体特性的描述非常细致，给观测者带来了很大的便利，但是对于遥远星系的描述这些软件还不是很完善。随着大型观测设备的出现，人们对未知天空的观测变得越来越遥远，星系模拟也成为评估大型天文设备很重要的一部分。文[1]介绍了SkyMaker，该软件很细致地描述了点扩散函数，包括大气扰动、望远镜运动、像差、衍射等各种因素，星系形态特征利用了De Vaucouleurs方法进行描述。文[2]利用Shapelet的方法对星系的形态进行模拟，该方法对星系形态采用多维模拟，文中与小波(Wavelet)等方法进行比较，对残差进行分析，证明其方法的可靠性。LSST[3]成立了专门的图像模拟团队，该团队主要工作是为LSST做模拟。他们的模拟针对全天区，利用宇宙模型建立星表确定天体的分布，并且建立了多种天体现象的模型(天体模型包括恒星、星系、小行星等)，在传输过程中考虑天气情况、大气扰动等，最后还考虑了结合巡天设备影响产生的点扩散函数。利用准确的模拟结果可以对LSST系统的各项性能进行分析处理。另外，LSST对天体形态的模拟采用蒙特卡罗的方法，模拟光子的射入。文[4]提到了模拟的两种不同方法，一种是类似LSST模拟光子射入的方法，另一种是按照像素模拟的方法，并比较两种方法的差异，文中用了简单的数学模型，主要强调处理的速度。

本文的模拟利用哈勃(Hubble)望远镜的数据和观测条件。图像模拟工作包含：(1)模拟天体的分布和天体的形态，主要是星系的形态；(2)模拟观测条件以及各种噪声，包括天光背景、暗流、读出噪声等；(3)模拟仪器的点扩散函数；(4)对观测结果进行统计分析。图1为图像模拟的流程图。

本文的主要目标是利用简单通用的模型对天文图像实现快速的模拟。在程序编写上应用更多的优化，使程序达到快速可用的目的。本文还利用SExtractor对模拟结果进行提取，将提取的结果与参考结果进行比较分析，通过分析确定模拟结果在实际应用上的可行性。

1 图像模拟的关键模型及方法

本文模拟利用HUDF ACS WFC i(F775)波段的数据，该数据生成的图像约为10 000×10 000像素，由于图像在空间上存在一定的旋转，图像4个角包含未拍摄的区域，为了实验方便，截取了其中4 500×4 500像素的中间区域，面积约为5.1 arcmin2。采用天文学中一些简单高效的模型，主要为了对数据进行快速处理。在观测条件上，模拟了哈勃望远镜在600 km高空的观测条件。

1.1 星系的模型

星系的模型采用Sersic模型[5]，数学表达式如下：

其中，I(R)为R处的流量密度；Ie为有效半径处的流量密度；Re为有效半径；n为Sersic系数，通过该系数可以控制星系的弥散程度；bn为常数。通过变换可以将(1)式变换成如下的表达式：

其中，I0为中心位置的流量密度；n为Sersic系数；a为常数，表示当R为a时的流量密度为中心位置流量密度的e-1倍。

图1 图像模拟流程图Fig.1 Flow Diagram of Image Simulation

若星系为圆形区域，对(2)式在圆内积分，得到星系的流量：

其中，Γ为gamma函数。

在模拟中，通常用椭率表示星系的形状，这就需要将星系的流量密度分布在整个椭圆平面内，定义椭圆的椭率为e，则(2)式中R可以用归一化的椭圆半径表示：

其中，x，y表示坐标，坐标的原点为椭圆的中心。

星系在图像中的表现，除了弥散程度和形状的特性外，还存在方向性，这就需要对已经计算出来的椭圆形的Sersic分布进行坐标变换，假设星系坐标系(长轴为x轴，短轴为y轴)与图像坐标系的夹角为θ，则变换矩阵为

新的坐标位置(x′， y′)为

通过(1)～(6)式的计算可以模拟星系的形状、大小、方向等要素，而对星系的模拟需要的参数，如n(Sersic系数)、Re(有效半径)、亮度等都是通过对哈勃超深空场ACS WFC的数据进行提取获得。在文[6]中利用SExtractor对哈勃超深空场的数据进行提取，并和文[7]提供的数据进行对比，结果比较理想，文中用到的参数均来自于文[6]提取的星表。

1.2 观测条件及噪声的模拟

对于观测条件的模拟，在图像中的表现其实就是对噪声的还原，还原的真实与否是影响模拟图像效果的最主要因素。观测产生的噪声主要包括固有噪声、信号噪声等。

固有噪声是指源于探测器自身的噪声，包括热噪声、散粒噪声、生成-复合噪声和闪变噪声，在文中用读出噪声表示。为了和哈勃的数据进行对比，读出噪声利用哈勃的曝光时间计算器[8]计算得出。

信号噪声主要包含背景噪声和由于信号的量子性产生的随机噪声。背景噪声利用了哈勃望远镜测得的天光背景噪声，包括黄道光和赤道光。

对于信号，由于其量子性，到达探测器会有一定的随机性，基于这种特性也会产生一定的噪声，在文中用泊松分布表示。

1.3 点扩散函数模拟

对于地面望远镜，影响点扩散函数的主要因素有大气扰动和仪器自身的光学特性。而对于哈勃望远镜这样的空间设备，观测条件要远远优于地面望远镜，因此对点扩散函数的影响仅仅来自观测设备的光学特性。

文中的模拟是对空间望远镜成像模拟，点扩散函数采用高斯(Gaussian)函数表示：

其中，I0为常数，表示中心位置的流量密度，可以控制高斯函数的幅度；[x，y]为图像坐标；σ用来控制图像的质量。对于高斯函数通常可以通过所给的半高全宽(FWHM)求得：

对于地面望远镜通常需要考虑大气扰动的效应，所以模拟中用到的点扩散函数需要考虑观测设备和大气扰动的综合效应，对于这类点扩散函数模拟通常采用Moffat[9]模型：

其中，I0为中心位置处的流量密度；α和β为尺度因子，可以控制函数的弥散程度；[x，y]为图像坐标。α和β关系如下：

文[10]中提到，当β＝4.765时，Moffat模型可以很好地模拟大气扰动。

1.4 子像素级处理

模拟的目的之一是对未来某一个观测设备的评估，文中模拟利用的数据只符合Hubble ACS的特性，例如像素的大小为0.03×0.03 arcsec2，为了使程序能够为其他观测设备评估利用，模拟需要做到子像素级。子像素级处理主要应用在天体形态、亮度的模拟过程中，根据需要将原像素分成若干个子像素，根据天体的特性，在子像素中计算得到该像素的值，计算量十分庞大，这就需要对程序进行优化。

1.5 并行计算

在计算天体的形态、亮度等要素时，由于将原有像素拆分为子像素计算，这样导致计算量倍增。另外，每一个天体都是独立的个体，单独计算对其他天体不会产生任何影响。所以，文中将计算每一个天体的形态、亮度等要素进行独立计算，这样，将模拟中的所有天体按照中央处理器的数目进行平均分配，不同的中央处理器计算的结果返给主进程进行合并，形成最终模拟的整幅图像。

2 结果分析

2.1 参考样本分析

利用文[6]提供的数据，其中包含了提取的Sersic系数，而Beckwith提供了一份更加准确的数据，不过在该数据中对于天体弥散程度的描述只给出了高斯参数，但是其位置和亮度更加准确。事实上，文[6]的数据也给出了一些与Bechwith数据的对比，例如，前者的数据给出了中心位置与Bechwith的对比，本文对两个数据做了更详细的对比。

图2是对两份数据提取天体亮度的对比，蓝色线为文[6]数据天体亮度的统计，红色线为Bechwith的数据天体亮度的统计。从图2可以看出，两份数据的亮度存在一定偏差，文[6]的数据天体亮度要比实际偏大。

图3是对两个样本有效半径的对比结果，横坐标是提取天体的亮度信息，纵坐标是Benchwich的数据中天体的有效半径与文[6]的数据中相应的天体有效半径之差。从图3可以看出，差值绝大部分集中在0值附近，但是文[6]的数据有部分半径偏大。

图2 Coe D.等提供的数据与Bechwith的数据天体亮度对比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]的数据模拟的图像与真实图像存在亮度偏差，另外，对于天体的大小也有部分天体与真实图像存在一定差别。从这一点，也可以解释后面章节中实验结果图4的模拟结果存在一些比原图大的天体的现象。

2.2 模拟结果及分析

图4展示了模拟的结果，图4(a)是模拟图像，图4(b)是经过MultiDrizzle程序处理后的ACS F775波段的图像。由于原图在空间存在一定角度，并且4个角上没有目标信息，所以本文用到的数据和图像是对原图的截取，截取了原图x方向3001到7500像素，y方向3001到7500像素的方形区域，总共4 500×4 500个像素，经过MultiDrizzle程序处理后的ACS F775波段的图像每个像素的视场约为0.03×0.03 arcsec2，模拟区域的视场为5.1 arcmin2。哈勃望远镜团队给出ACS WFC i(F775)波段零点星等(zeropoint)的值为25.654，这里零点的意思是流量为1e-时天体的亮度值。根据(11)式可以计算天体的流量：

由(11)式计算的流量值是信号的流量。而噪声包括天光背景、暗流、读出噪声等，是利用哈勃望远镜的曝光时间计算器计算得到。计算结果以e-为单位的流量。所有计算的结果除以增益(单位，e-/DN)得到表现在图像上的数字值。

图4 模拟图像结果Fig.4 Simulation Image vs.Original Image(a： Simulation Image， b： Original Image)

HUDF ACS WFC I波段的极限星等为29等，所以本文提取的截止亮度到29等。文中用SExtractor作为提取工具，SExtractor的提取主要参数见附录。表1是对提取数目的统计，参考的对象为文[6]的数据。表中有4项数据，第1项是对原图提取数据的统计，第2项是对模拟图像提取数目的统计，第3项是粗略地做了匹配统计，匹配的方法是将模拟数据提取的天体中心与文[6]数据的天体中心作比较，选出距离最接近的一对被认为是同一个天体。第4项是计算匹配率的比例，该比例通过(12)式计算得到：

其中，rmatch为匹配率，Nref为参考数据提取天体的数目；Nmatch为对模拟数据提取天体的数目。从表1匹配的数目可以看出，这种针对位置的匹配度均在85%以上，考虑到噪声影响和提取方法的差异，这种匹配程度是可以接受的。

表1 利用模拟图像提取的天体与HUDF catalog天体匹配数目统计Table 1 Extracted galaxies analysis

图5是对小于29等的天体数目的统计，从图5可以看出，对于模拟图像提取的天体数目和参考数据天体的数目随亮度变化的趋势是比较一致的。

图6是将模拟数据提取的天体有效半径与参考数据天体的有效半径进行对比，采用了计算匹配天体有效半径之差的方法进行对比，匹配天体的方法见表1，在这里只是确定了位置上的匹配，这两个统计的样本都是亮度小于29等的天体。图6(a)统计该差值随着亮度分布的变化，横坐标是亮度，纵坐标是参考数据中天体的有效半径与相应的模拟数据提取的天体有效半径之差，从分布图上可以看出，亮度越小，有效半径的差值越小。图6(b)是对差值的分布进行统计，从统计可以看出，这一差值绝大多数分布在0值左右，有部分存在较大的偏差。导致这种差异的原因有很多，其中包括提取样本的精度(2.1节中对两种参考样本的分析)，图像模拟中加入的信号噪声、背景噪声和设备噪声及利用SExtractor对天体进行提取的参数也存在一定差异。

图7是对模拟图像提取的天体亮度与参考数据的天体亮度进行对比，对比的方法同图6(a)，通过亮度之差随天体亮度的分布进行分析。图中，横坐标是亮度，纵坐标是参考数据中天体的亮度与相应的模拟数据中提取的天体亮度之差。统计的样本依然是亮度小于29等天体。从图7可以看出，这一亮度差值大部分分布在0轴上下，偏差绝大多数集中在0.5以下。从图7还可以看出一个规律，对于亮度越大的天体，模拟后再提取，其亮度值越准确。

图5 提取天体数目统计对比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 性能分析

本文中图像模拟软件以及统计的数据生成都是用JAVA程序开发。在天文中用JAVA开发的程序较少，本文之所以选择JAVA程序开发主要因为JAVA程序存在以下优势：(1)JAVA去掉了指针，有自动的垃圾回收机制，这样减小了编写程序时由于内存泄露出现的问题，因此，能够将更多的时间用在算法的编写上；(2)JAVA程序是面向对象的编程，在程序的可读性上要好于C程序。(3)JAVA支持并行计算，且编写简单，对于本程序的并行计算有很大的帮助。

另外，本文的程序运行在Linux系统下，硬件配置CPU为Intel core i5 4核3.1 GHz，内存为16 GB，在程序运行时，4核全部使用，计算天体的数目为4 107，计算所需要的时间为1 267 s，平均每个天体所需要时间为308 ms。程序运行时峰值所需要的内存为4.1 GB。

程序运行的大部分时间和内存都是消耗在利用Sersic模型计算天体的形态，在计算天体的形态时，由于1.4节所述子像素的需求以及计算的精度，需要将一个像素分解成若干个小像素，计算时间和内存的消耗与像素分解的个数成正比。

对于星系形态的计算采用了二维的数值积分，为了计算更加准确，本文将一个像元分解为若干个子像元，通过利用子像元积分得到每一个像元的准确数据，文中兼顾性能需求没有对像元进行无限划分，只是分解为32×32个子像素，4 500×4 500像素的图像就变成了144 000×144 000像素的图像，这样数据量已经很大，采用并行及算法上的优化使得每个天体的计算时间达到308 ms。

3 总结及展望

本文对图像模拟用到的方法都是天文中比较基础的模型，例如Sersic模型、高斯模型。将这些基础模型应用到模拟中，可以达到简单高效的目的。但是简单高效只是在模拟结果真实正确的基础上才能成立。通过本文第2部分的分析可以发现，模拟的结果与真实数据存在一定差异，但是基本特征都能够正确地表现(例如大小、方向、亮度等)；通过统计分析可以看出，对于模拟图像提取的天体的数目、亮度、大小与参考数据比较，错误的数量都在可以控制的范围内。另外，本文对模拟程序实现了并行计算，而且在程序中已做了一些优化(例如，三角函数用插值方法计算)，通过对性能的分析可以看到，计算天体的平均时间达到了毫秒级，再次说明了本文的模拟成像速度快。

模拟程序中，点扩散函数只用到高斯函数，虽然简单，但是对于反映图像的真实性还存在一定的差别。另外，模拟的目的主要是为一些新的巡天项目分析服务，所以下阶段的目标是在模拟程序上对点扩散函数进行改进，加入真实的光学系统特性。另外需要引入光谱的模拟以适应多波段成像模拟。同时，在分析上，将模拟与新的巡天项目相结合，得到与不同硬件相关的分析。