射电天文 / 观测中的一些基础知识(4)：简单图像处理与流量换算

我是一个梦蛋

去年讲了一些射电天文的基础知识
然后挖了个坑讲厘米波谱线数据处理举例
这次开始介绍一些新内容。
推荐阅读
天文中的常用Python操作（FITS） (astrogeeker.com)  Charon学姐的网站
五、对射电cube及图像的常见简单处理

我在gitlab上放了个包，是我大三时候写的简单射电图像包，比较初级，欢迎来魔改源码xs
astroR2 / hiviewer · GitLab
注：本文中的草履虫(划掉)M33星系为FAST早期测试数据(的一部分，可能有点问题)，仅供教学科普使用。
1. 基本信息

我们前面提到了一个data cube是三维的，二维空间+一维速度(频率)，数组的每一个点表示强度。下面我从hiviewer的示例截图说明。
天文上常用的数据格式是fits，我们可以使用FITS File Handling (astropy.io.fits) — Astropy v5.2的方法读取。
读取header：


数组维度和wcs。World Coordinate System (astropy.wcs) — Astropy v5.2
注意维度是(z, y, x)，或者说(velo, dec, ra)


使用Spectral Cube documentation — spectral-cube v0.6.0查看cube


速度带astropy.units的单位 Units and Quantities (astropy.units) — Astropy v5.2


波束用这个包定义Radio Beam — radio-beam v0.3.4


坐标转换主要根据SkyCoord — Astropy v5.2和WCS的all_pix2world、all_world2pix函数，因为像素坐标和天球坐标的转换还是比较常用的。


已经有很多成熟的Python包可以用了2333
2. 切一片看看

三维cube固定一个速度切一刀，我们便可以看到空间截面。



正宗俄罗斯大列巴，来自网络

啊，不是这个，是切星系



M33的一个截面，RA和Dec是赤经赤纬

如果固定RA或者Dec切一刀，那么得到了一种pv图，位置position-速度velocity图。这样我们能看到HI盘在不同位置处的速度。



固定ra切，横着的一长条是银河系，速度在0km/s附近


3. moment maps

直译矩图，但是感觉怪怪的。它们在射电中非常非常常用，定义是
M_0= ∫I_v  \rm{d}\it{v}
M_1= \frac{1} {M_0} ∫v I_v  \rm{d}\it{v}         
M_2= \frac{1}{M_0}\int (v-M_1 )^2 I_v  \rm{d}\it{v}
其中 I_v 代表辐射强度， v 代表速度.
M_0 是零阶矩，流量密度对速度的积分，其含义就是我们常说的流量图。
M_1 是一阶矩，可以看出来它就是速度按辐射强度加权的平均(联想到了热统的分子运动平均速度)，自然是速度的量纲。速度图告诉我们每条谱线按辐射强度加权后的中心速度。
M_2 是二阶矩， (v-M_1 )^2 代表速度的弥散，也是按强度加权的平均，速度平方的量纲。于是可以定义N阶矩。对moment2开根号得到sigma map，乘个系数 2\sqrt{2\ln2} 就是FWHM图，告诉我们每条谱线加权后的线宽。
例如下图是CO的一篇文章中的M0,M1,M2
Moment maps of the C18O J = 3–2 line. Left: Moment 0 map, integrated... | Download Scientific Diagram (researchgate.net)


我们可以使用Moment maps and statistics — spectral-cube v0.6.0中的函数方便的截取你感兴趣的速度区间做moment maps。

4. p-v图

当然我们也可以在RA-Dec平面上任意取点切一刀，这在射电上研究气体stream，有没有联系等等很常用。
以沿主轴方向的pv图为例。首先对星系Moment0的等高线拟合一下，找到中心和主轴。这里我使用的方法是下面第10条



红色直线就是主轴了

当然测光里也有现成的方法Elliptical Isophote Analysis (photutils.isophote) — photutils 1.6.0
知道椭圆的中心和主轴之后，我们便可以愉快地切pv图了。这里R2参考了多波段的任老师代码魔改了一下，感谢~



有旋转曲线的feel了

5. 抽一条谱线看看

好拔，才想起介绍这个，不过很简单，抽取某个像素坐标里的谱线。记得cube的shape是(z, y, x)，这样从data[:, y, x]取谱线



多峰谱线。左边是M33，右边是银河系

6. 卷积和重新投影

现在我们有FAST, Arecibo和VLA的三个cube，想比较一下它们的效果，但是我们注意到它们的分辨率和像素大小不一样，这时便需要统一一下。
我们可以使用Smoothing — spectral-cube v0.6.0里的convolve_to函数将它们卷积到最差的分辨率上(这里是Arecibo 3.5角分)，使用重新投影的方法将它们reproject到最大的格点上(这里是FAST，1.5角分)，这样便可以逐一像素的比较了。重新投影的方法可以参考以下三种reproject方法

• 一个mosaic软件Montage - Image Mosaic Software for Astronomers (caltech.edu)和它的python wrapper包 Python Montage wrapper (montage_wrapper) — montage_wrapper v0.9.9.dev237 (montage-wrapper.readthedocs.io)
  
• 荷兰人的天文图像包Kapteyn Package for Python 3 — Kapteyn Package (home) (rug.nl)
  
• 以及一个就叫reproject的包Image reprojection (resampling) — reproject v0.9.1
  

当然如果觉得麻烦，欢迎使用我打包好的函数(虽然大三时候的R2就是调包侠)
例如下图是的草履虫，灰色线是模糊后的VLA，绿色是Arecibo的等高线。


更高级的好像在知乎写不下，不过这些常用的初级射电操作差不多了。
六、其它可能用到的

1. 流量的单位换算

这是我想写这个系列的初衷，但是到这里第四篇才想起来写[呆]


我们在第一篇讲到了流量密度的定义
S_\nu = \int_{\Omega_S}I_\nu(\theta,\phi)\cos \theta d\Omega
I_\nu是辐射强度。积分是对立体角的的积分。I_\nu的单位就是 \rm{Jy·sr^{-1}} ，sr是球面度。
观测的时候是我们认为直接测量的是望远镜的一个波束(beam)收到了多少Jansky的流量，所以这个流量 S 与波束大小有关。所以我们关心的是源的辐射强度，这是强度量，应该与用什么波束大小的望远镜无关的。
假设波束立体角为 \Omega_b ，源的强度为 I ，流量简单的写为 S=I\Omega_b 来看看量纲。
单位就是 \rm{Jy} = \rm{Jy·beam^{-1}} \times \rm{beam} ，I 的单位就是 \rm{Jy·beam^{-1}}
这挺合理，Jy per beam是文献中很常见的辐射强度单位。

我们的观测数据如果已经做成了cube，里面每个像素点又是什么单位呢？
这时看一看fits的header，BUNIT里如果是Jy/beam，K都常见(温度也可以转化嘛)
你可能想把源的流量在空间平面上加起来算个和，但是这时要注意单位转换：
\rm{Jy·beam^{-1}} \rightarrow  \rm{Jy·pixel^{-1}}
举个栗子，beam的FWHM是3角分，pixel是1角分的情况。怎么转换呢？



思考.jpg

首先我们知道一个理想的二维椭圆高斯的立体角怎么求:
一维情况大家都知道 \begin{equation}\int_{-\infty}^{\infty} \exp{-\frac{x^2}{2\sigma_x^2}} dx= \sqrt{2\pi}\sigma_x\end{equation}  
由Gaussian Function -- from Wolfram MathWorld
x_{0}= \sigma_x \sqrt{2 \ln 2} (这个很好算，就峰值一半对应的x)
所以 \begin{equation}FWHM=\theta_x=2\sqrt{2\ln2}\sigma_x\end{equation}  
二维情况，xy独立，这个对全空间的面积分就是我们要求的立体角
\begin{align}I_{all}&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \exp{[-\frac{x^2}{2\sigma_x^2}-\frac{y^2}{2\sigma_y^2}]} dxdy \\ &= 2\pi\sigma_x\sigma_y \\ &= \frac{\pi\theta_x\theta_y}{4\ln2}=\Omega\end{align}  
这里利用了 \sigma_x\sigma_y=\frac{\theta_x\theta_y}{8\ln2}

好了，波束立体角现在是 \Omega_b = \frac{\pi\theta_x\theta_y}{4\ln2} 。要转换这两个单位，只差一个波束立体角，我先放结果：
S(\rm{Jy·pixel^{-1}})= \it{I}\rm{(Jy·beam^{-1}}) / \Omega_b’
这里的 \Omega_b' 是取像素单位的。对于FAST，我们取 \theta_x，\theta_y 都是 3' 先。对于pixel是1角分的cube，FWHM可以用3个pixel来代入，这样得到的立体角大概是10，就是一个beam对应10pixel大小。
下面来解释一下：


假设源的辐射强度很均匀，就像左边的红色数值都是10 \rm{mJy·beam^{-1}}(我们记得它是强度量，就像压强一样，单位面积的压力).
而 \rm{mJy·pixel^{-1}} 是可以求和的，你可以对很多pixels求和算一个总的flux(我们记得它是广延量，就像压力一样，你可以求和)
所以我们只需要一个波束立体角对应多少个pixel的转换系数， \Omega_b’
10 \rm{mJy·beam^{-1}}除以一个beam包含的pixel数(10)之后变成右边的橙色1 \rm{mJy·pixel^{-1}}(当然我只画了9个方格，那个面积分是10，明白意思就行)。
这样10个橙色的1mJy加起来便还原了原来一个beam中的10mJy。或者1个pixel里的流量1mJy是原来一个beam中的10mJy的十分之一。

现在就可以对 \rm{Jy·pixel^{-1}} 单位的流量在空间平面求和了
S_{\mathrm{sum}}(\rm{Jy})=\sum_{source}^{N}{S(\rm{Jy·pixel^{-1}})}
当然如果还有谱线方向，类似于求moment0
S_{\mathrm{sum}}(\rm{Jy·km\ s^{-1}})=\sum_{source}^{N}{\it{S}\rm{(Jy·pixel^{-1}})} \Delta \it{v} \rm{(km\ s^{-1})}
或者参考[1]的写法
S_{\mathrm{sum}}=\frac{\Delta z}{\Omega_{\mathrm{PSF}}} \sum_i S_i
这样就包含了单位换算。 \Omega_{\mathrm{PSF}} 是就是上面说的像素单位的beam solid angle。

这个转换当时困扰了我一会儿，有些其它的参考可能也有用：
科学网—Jy/beam到Jy/pixel的转化 - 张国印的博文 (sciencenet.cn)
radio astronomy - Converting Jy/beam to Jy? - Astronomy Stack Exchange
convolution and regridding 一个之前觉得很好的blog但是打不开了
2. 我还不知道写啥

放一下之前零零碎碎的吧，比如射电软件
以及没啥人看的射电科普回答
R2学射电也没几年，文章有错误欢迎指正。感谢
估计射电天文基础系列就暂时告一段落了，直到下次我知道写点什么水一篇。寒假期间可能蹭蹭三体的热度，回顾一下地球往事三部曲中出现了哪些天文学家，毕竟是刘电工的小说让我觉得天文学家可以毁灭/拯救世界！
啊，是遥远的高中回忆，也可以水一篇。
<hr/>封面图：  Arecibo Message  ，唉，缅怀
欢迎关注R2的专栏：予日行辰，不定时更新。。。
本文作者：R2
本文审核：不愿透露姓名的喵同学
本文所使用图片审核：不愿透露姓名的师兄


分类整理：天文、Linux、编程、杂谈、科幻、科研
参考

• ^SOFIA2 Manualhttps://github.com/SoFiA-Admin/SoFiA-2/wiki/documents/SoFiA-2_User_Manual.pdf