Chandra图像定标

伪彩色这是2008年Paolo Tozzi 来北京做的X射线数据处理讲座的笔记,一直拖到现在才开始整理。这里是第一部分,主要是些准备工作。高能数据处理要用HEASOFT,Chandra数据还要用到Ciao和CALDB,按照说明文件一步步安装并不困难,只是2G多文件下载起来要费些时间。附加脚本也经常用到,解压到安装路径下就可以了。安装好之后,先用 heainit 和ciao初始化运行环境,若没有错误提示就可以运行下面的命令了。

如果手头上没有数据。可以先到Chandra数据中心找个自己喜欢的目标,数据通常会在保密期过后自动公开,既避免了对同一目标的反复观测,也能最大限度地发挥已有数据的价值。 我这里选择了红移0.3的 MACSJ0257.6-2209 (即Abell 402, Obs ID 3267 )。输入到左上角的Target Name中,点Search,会得Chandra已有的观测数据列表。我们选个曝光时间短点的,数据文件小些,加到获取列表中(Add to Retrieval List)。点击retrieval之后,服务器会自动将文件打包放到匿名ftp的临时目录下。也有一个基于python的自动下载脚本。 官方详尽的下载介绍由此跳转

压缩包里东西很多,图像处理只需要部分文件:
primary目录下,记录坏点的 bpix,确定指向的 asol
secondary目录下,观测事例 evt ,描述可用时间间隔的 flt,储存观测参数的 pbk 文件,
后面的数字代表处理级别,0为原始数据,1为基本处理数据,2为科学数据
完整命名规则由此去

为防止误操作原始数据,简化文件名,先做符号链接

gunzip *.gz
ln -s *evt1.fits evt1.fits
ln -s *bpix1.fits bpix1.fits
ln -s *asol1.fits asol1.fits
ln -s *flt1.fits flt1.fits
ln -s *pbk0.fits pbk0.fits
ls -l

启动ciao,用dmlist检查fits文件头中的基本信息。检查文件区块

$dmlist  evt1.fits blocks

--------------------------------------------------------------------------------
Dataset: evt1.fits
--------------------------------------------------------------------------------

     Block Name                          Type         Dimensions
--------------------------------------------------------------------------------
Block    1: PRIMARY                        Null
Block    2: EVENTS                         Table        16 cols x 597013   rows
Block    3: GTI3                           Table         2 cols x 2        rows
Block    4: GTI2                           Table         2 cols x 3        rows
Block    5: GTI1                           Table         2 cols x 2        rows
Block    6: GTI0                           Table         2 cols x 2        rows
Block    7: GTI6                           Table         2 cols x 2        rows

可以看到 fits 文件实际是按表格组织的,GTI是 可用时间间隔(Good Time Interval),记录起止时间每次两列。光子事件共有 597013 行,16列的具体含义可以用”dmlist evt1.fits cols”查看。
在ACIS的10片ccd中有两片是后照式BI(S1和S3),处理方式略有不同。详细介绍参考申请指南中的相关页面 下面以ACIS-I的四片FI式CCD为例简要介绍观测数据的定标处理过程。

dmcopy "evt1.fits[ccd_id=0,1,2,3]" evt1_acis_I.fits

将CCD数据提取到新文件中,然后用acis_process_events作初步的处理。

$ dmlist evt1.fits header | grep DATAMODE  | awk '{print $3}'
VFAINT  # 检查观测模式, 确定下面的 check_vf_pha 参数
$ dmlist evt1.fits header | grep Focal | grep FP_TEMP| awk '{print $3}'
153.44601440  # 检查焦面温度,153K需要做 CTI 修正, 对应下面的 apply_cti 参数

下面用 punlearn 更新命令参数,虽然也可以直接在命令后加参数,不过前者看起来比较清楚。

punlearn acis_process_events
pset acis_process_events infile=evt1_acis_I.fits
pset acis_process_events outfile=evt1_updated.fits
pset acis_process_events eventdef=")stdlev1"
pset acis_process_events check_vf_pha=yes
pset acis_process_events apply_cti=yes
pset acis_process_events rand_pha=yes
pset acis_process_events doevtgrade=yes
pset acis_process_events apply_tgain=yes
pset acis_process_events gradefile=CALDB
pset acis_process_events gainfile=CALDB
pset acis_process_events calculate_pi=yes
pset acis_process_events badpixfile=bpix1.fits
pset acis_process_events acaofffile=asol1.fits
pset acis_process_events stop=none
pset acis_process_events mode=h
pset acis_process_events verbose=5
pset acis_process_events clobber = yes
acis_process_events

Chandra的正常光子事件是在可用时间间隔内 grade 为 [02346]之一,状态为0的那部分。我们用dmcopy来提取。

dmcopy "evt1_updated.fits[EVENTS][grade=0,2,3,4,6][status=0][@flt1.fits]" evt2.fits

这就是从evt1到evt2的标准处理流程,如果没有什么特殊要求的话,直接用primary目录下的evt2文件也是一样的。在ds9中检查一下,应该和压缩包中的 full_img 图像是一样的。

为避免穿越辐射带,遭遇太阳风暴等情况带来的高能背景涨落,还要选取背景信号稳定的观测时段。下面统一处理四个CCD,谨慎一点的话可以对他们分别处理。
先在ds9中手动选出包括星系团在内的所有亮源,保存成souce.reg文件。然后扣除这些区域提取背景光谱,以免X射线源流量变化带来的干扰。接下来观察背景光子计数率变化曲线,时间间隔取100秒。

punlearn dmextract
pset dmextract "evt2.fits[exclude sky=region(source.reg)][bin time=::100]"
pset dmextract outfile=check.pha
pset dmextract clobber=yes
pset dmextract opt=ltc1
pset dmextract mode=h
pset dmextract verbose=5
dmextract

借助现成的 S-Lang 脚本,详细步骤可参考此页面

slsh
slsh> require ("lightcurves");
slsh> lc_sigma_clip("check.pha");
....

GTI limits calculated using a count-rate filter:
  (count_rate>1.2454 && count_rate<2.01135)

The corresponding times are:
  ((time >= 121965879.253540) && (time < 121986779.253540)) ; 20.47 ksec, bin 1

....
slsh>quit

程序默认剔除所有计数率偏离均值3σ以上的异常时段。

观察背景光变

我们把脚本输出的正常时段保存成新的gti文件来更新源文件,注意要删除其中的等号,否则报错。

punlearn dmgti
pset dmgti infile=evt2.fits
pset dmgti outfile=evt2.gti
pset dmgti userlimit = "  ((time > 121965879.253540) && (time < 121986779.253540))"
pset dmgti clobber=yes
pset dmgti mode=h
pset dmgti smooth=yes
dmgti

更新fits文件

dmcopy "evt2.fits[@evt2.gti]" evt2_clean.fits

不放心的话可以再重新提取背景光谱检查一下。如果是对几片CCD分别处理,还要对其余三片重复上述过程,然后再把它们重新合并起来

punlearn dmmerge
pset dmmerge infile=ccd0_clean.fits, ccd1_clean.fits, ccd2_clean.fits, ccd3_clean.fits
pset dmmerge outfile=evt2_clean.fits
pset dmmerge outBlock=events
pset dmmerge mode=h
pset  dmmerge clobber=yes
dmmerge

如果有多次观测可以在分别处理后通过merge_all合并。定标就算完成了,接下来终于可以提取数据了。在命令dmcopy中,可以通过 energy 来选取波段,用bin来选取或合并像素

dmcopy "evt2_clean.fits[events][energy=500:5000][bin x=2000:5500:2,y=2000:5500:2]" img_tot.fits

为了更好地展示图像,可将不同的能段赋予不同颜色,产生伪彩色图。这里我们选取0.3-1 keV 为红色,1-2keV 为绿色,2-5 keV 为蓝色

dmcopy "evt2_clean.fits[energy=300:1000][bin x=2000:5500:2,y=2000:5500:2]" img_cl_red.fits
dmcopy "evt2_clean.fits[energy=1000:2000][bin x=2000:5500:2,y=2000:5500:2]" img_cl_green.fits
dmcopy "evt2_clean.fits[energy=2000:5000][bin x=2000:5500:2,y=2000:5500:2]" img_cl_blue.fits

在Chandra的高角分辨率下X射线光子分布十分弥散,我们可以用aconvolve命令做一下高斯平滑,方便查看。

punlearn aconvolve
pset aconvolve clobber=yes
pset aconvolve infile=img_cl_green.fits
pset aconvolve outfile=img_cl_greensm.fits
pset aconvolve kernelspec="lib:gaus(2,5,1,2.5,2.5)"
pset aconvolve method=fft
pset aconvolve verbose=5
pset aconvolve mode=h
pset aconvolve normkernel = area
aconvolve

上面这步对三个波段都要做,我这里就不重复了。接下来用dmimg2jpg合并三色图像就可以得到伪彩色图了。

punlearn dmimg2jpg
pset dmimg2jpg infile=img_cl_redsm.fits
pset dmimg2jpg greenfile=img_cl_greensm.fits
pset dmimg2jpg bluefile=img_cl_bluesm.fits
pset dmimg2jpg outfile=color.jpg
pset dmimg2jpg maxred=0.7 maxblue=0.7 maxgreen=1.5
pset dmimg2jpg minred=0 minblue=0 mingreen=0
pset dmimg2jpg gridsize=60
pset dmimg2jpg fontsize=1
pset dmimg2jpg showaimpoint=no
pset dmimg2jpg clobber=yes
pset dmimg2jpg scalefunction=power
pset dmimg2jpg scaleparam=0.6
pset dmimg2jpg mode=h
dmimg2jpg

伪彩色

当然,生成jpg格式主要是方便检查核对。图片确实不够漂亮,但包含了很多有用的信息。除了明亮的星系团,图中还有很多颜色各异的亮点,通常是AGN。蓝色的温度较高,由于低能段被介质吸收,因而缺少红色,不然便是白色。而红色的温度较低,高能光子较少。更具体的情况可以直接查看X射线光谱,将在下篇中介绍。

订阅评论
提醒
11 评论
最旧
最新 最多投票
内联反馈
查看所有评论
firefly
2010 年 8 月 7 日 18:45

你好,看了用《matlab做星点识别》文章后,很有启发。我现在刚接触星点识别这个工作,就是要用matlab编程实现星图中星点的定位,要定位到亚像素级。现在用fitsread函数读出图像的数据后,想用matlab显示出来,并且对图像进行滤波去噪处理,不知道该怎么做。还有,fitsread读出的图像矩阵的值是什么含义,能当作灰度图像的灰度值那样处理吗?还往能不吝赐教,非常感谢^_^

hbcai
2010 年 8 月 9 日 10:05

整理得好。期待继续。

胖麻雀
2010 年 12 月 18 日 09:25

为什么在你的处理流程里没有对afterglows的修正?

SN1987A
2011 年 6 月 21 日 14:13

先在ds9中手动选出包括星系团在内的所有亮源,保存成souce.reg文件。然后扣除这些区域提取背景光谱,以免X射线源流量变化带来的干扰。

pset dmextract “evt2.fits[exclude sky=region(source.reg)][bin time=::100]”

张鹏飞
2012 年 3 月 27 日 16:00

你好!我正在学x-ray的处理,不知道chandra,XMM-newton,RXTE的那个简单,我想先学个简单的!这三个卫星哪个好一点?你能给我个答复吗?

菲菲
回复给  gerry
2012 年 3 月 28 日 10:50

多谢了,以后有问题了,还需要你多多帮忙啊!