UFO中文网

 找回密码
 注册会员
查看: 151|回复: 0
收起左侧

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

[复制链接]
online_member 发表于 2023-1-14 18:56:18 | 显示全部楼层 |阅读模式
去年讲了一些射电天文的基础知识
然后挖了个坑讲厘米波谱线数据处理举例
这次开始介绍一些新内容。
推荐阅读
天文中的常用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:

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算560 / 作者:我是一个梦蛋 / 帖子ID:104513
数组维度和wcs。World Coordinate System (astropy.wcs) — Astropy v5.2
注意维度是(z, y, x),或者说(velo, dec, ra)

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算772 / 作者:我是一个梦蛋 / 帖子ID:104513
使用Spectral Cube documentation — spectral-cube v0.6.0查看cube

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算504 / 作者:我是一个梦蛋 / 帖子ID:104513
速度带astropy.units的单位 Units and Quantities (astropy.units) — Astropy v5.2

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算389 / 作者:我是一个梦蛋 / 帖子ID:104513
波束用这个包定义Radio Beam — radio-beam v0.3.4

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算950 / 作者:我是一个梦蛋 / 帖子ID:104513
坐标转换主要根据SkyCoord — Astropy v5.2和WCS的all_pix2world、all_world2pix函数,因为像素坐标和天球坐标的转换还是比较常用的。

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算54 / 作者:我是一个梦蛋 / 帖子ID:104513
已经有很多成熟的Python包可以用了2333
2. 切一片看看

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

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算169 / 作者:我是一个梦蛋 / 帖子ID:104513

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

啊,不是这个,是切星系

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算300 / 作者:我是一个梦蛋 / 帖子ID:104513

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

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

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算849 / 作者:我是一个梦蛋 / 帖子ID:104513

固定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)

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算969 / 作者:我是一个梦蛋 / 帖子ID:104513
我们可以使用Moment maps and statistics — spectral-cube v0.6.0中的函数方便的截取你感兴趣的速度区间做moment maps。

4. p-v图

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

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算624 / 作者:我是一个梦蛋 / 帖子ID:104513

红色直线就是主轴了

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

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算876 / 作者:我是一个梦蛋 / 帖子ID:104513

有旋转曲线的feel了

5. 抽一条谱线看看

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

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算204 / 作者:我是一个梦蛋 / 帖子ID:104513

多峰谱线。左边是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的等高线。

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算528 / 作者:我是一个梦蛋 / 帖子ID:104513
更高级的好像在知乎写不下,不过这些常用的初级射电操作差不多了。
六、其它可能用到的

1. 流量的单位换算

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

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算216 / 作者:我是一个梦蛋 / 帖子ID:104513
我们在第一篇讲到了流量密度的定义
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角分的情况。怎么转换呢?

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算421 / 作者:我是一个梦蛋 / 帖子ID:104513

思考.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大小。
下面来解释一下:

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算554 / 作者:我是一个梦蛋 / 帖子ID:104513
假设源的辐射强度很均匀,就像左边的红色数值都是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
本文审核:不愿透露姓名的喵同学
本文所使用图片审核:不愿透露姓名的师兄

射电天文 / 观测中的一些基础知识(4):简单图像处理与流量换算867 / 作者:我是一个梦蛋 / 帖子ID:104513
分类整理:天文、Linux、编程、杂谈、科幻、科研
参考


  • ^SOFIA2 Manualhttps://github.com/SoFiA-Admin/SoFiA-2/wiki/documents/SoFiA-2_User_Manual.pdf
您需要登录后才可以回帖 登录 | 注册会员

本版积分规则

UFO中文网

GMT+8, 2024-11-5 21:57

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表