天体运动的可视化及射电天文光谱数据处理

科学技术原理

Python 以及衍生的一系列数据处理分析软件包已经在天文领域获得了广泛使用,在天文数据处理和分析、数据管理、天文数据绘图、高性能计算、机器学习和深度学习、乃至望远镜管理等诸多领域,都有着广泛的使用和不俗的性能。这里着重介绍可视化的应用,即充分利用 matplotlib 库。本次任务会介绍相对简单的两种表现:第一,模拟天体运动轨道并可视化;第二,射电天文光谱数据可视化。

1. 天体运动可视化:通过物理公式计算模拟轨道的参数并将其显示在图像上

2. 射电天文光谱数据可视化:通过简单的数据提取和数据处理将结果显示在图像上1

设计方案

天体运动可视化

关于物理公式

A天体(质量为m)在与B天体(质量为M)作用下的万有引力的影响下按一定轨道运动。 如果我们用数学公式表示A天体此时接收到的力:

\vec{f} = -\frac{GMm}{|\vec{r}|^2}\frac{\vec{r}}{|\vec{r}|}

运动方程:

m\frac{d^2\vec{r}(t)}{dt^2}=-\frac{GMm}{|\vec{r}|^2}\frac{\vec{r}}{|\vec{r}|}

初始位置设置在x轴上,给一个y轴方向的初速度:

\vec{r}(t=0)=(\,R_0\,,\,0\,) \quad \quad \vec{v} (t=0)=(\,0\,,\,v_{y0}\,)

将运动方程改写成一阶微分方程(要使用使用四阶 Runge-Kutta 法求解):

\frac{d\vec{r}}{dt}=\vec{v} \quad \frac{d\vec{v}}{dt}=-\frac{GM}{|\vec{r}|^2}\frac{\vec{r}}{|\vec{r}|}

为了在程序上便于操作,将变量归纳为矩阵形式

X(t):=(\,\vec{r}\quad \vec{v}\,)=\left( \matrix{ x & v_x\\ y & v_y } \right)\\ F(X(t))=(\,\vec{v}\quad -GM·\vec{r}/|\vec{r}|^3)=\left( \matrix{ v_x & -GM·\vec{r}/|\vec{r}|^3\\ v_y & -GM·\vec{r}/|\vec{r}|^3 } \right)\\

因此把一阶微分方程改写成如下形式:

\frac{dX(t)}{dt}=F(X(t))

关于算法(四阶Runge-Kutta法)2

四阶 Runge-Kutta 法用于求常微分方程的较高精度的数值解。算法表示如下:

(Δ_1,Δ_2,Δ_3,Δ_4):=(0,Δt/2,Δt/2,Δt)\\ (β_1,β_2,β_3,β_4):=(1,2,2,1)\\ K_0:=O \quad \quad K_i:=F(X(t)+Δ_iK_{i−1})\\ X(t+Δt)=X(t)+\sum_{k=1}^4β_kK_k\\

其中O是零矩阵

初始值:

X(t=0)=\left( \matrix{ R_0 & 0\\ 0 & v_{y0} } \right)\\

射电天文光谱数据处理

由于实验数据格式的特殊性,需要用到额外的库进行数据处理

import astropy      # 天文数据处理
import spectral_cube  # 高光谱模块

并且数据格式固定,每一条数据都能单独被提取,因此实际操作并不复杂,只需要将数据对应到图像上即可

创新性描述

通过 Python 的可视化库 matplotlib ,不仅可以将天体辐射数据模型以一种更直观的方式展现出来,而且可以通过纯数学方式模拟展示一些天体运动的模型。只通过几行代码,就可以生成图表、直方图、功率谱、条形图、误差图、散点图等。 因此 matplotlib 库函数的运用在计算机视觉等专业上有着非同寻常的意义。

运行方法和参数设置

天体运动可视化

四阶 Runge-Kutta 法的代码表示:

while t < tmax:
 sum_x = np.zeros((2, 2))
 for k in [1, 2, 3, 4]:
 k_rk[k, :, :] = f(G_pa, M_pa, x_mat + delta[k]*k_rk[k-1, :, :]) # Ki计算
 sum_x += beta[k]*k_rk[k, :, :] 
 x_mat += (dt/6)*sum_x
 if count % cat_val == 0:
 x_lis.append(x_mat.copy())
 Nt += 1
 count += 1
 t += dt

呈现图像所需要的参数:

# 位置
 r_val = (x_t[i]**2 + y_t[i]**2)**0.5

 # 运动径向的单位向量
 er = np.array([x_t[i], y_t[i]])/r_val

 # 万有引力
 m_pa = 1.0
 f_vec = -1.0*(G_pa*M_pa*m_pa/(r_val**2)) * er.copy()
 plt.quiver(x_t[i], y_t[i], f_vec[0], f_vec[1], scale=vec_size**-1, color='red')

 # 速度
 plt.quiver(x_t[i], y_t[i], vx_t[i], vy_t[i], scale=vec_size**-1, color='blue')

 # 速度大小
 v_val = (vx_t[i]**2 + vy_t[i]**2)**0.5

 # 离心力
 m_pa = 1.0
 f_vec = (m_pa*v_val**2)/r_val * er.copy()
 plt.quiver(x_t[i], y_t[i], f_vec[0], f_vec[1], scale=vec_size**-1, color='green')

 # 轨迹
 plt.plot(x_t[0:i], y_t[0:i], '--')
 plt.scatter(x_t[i], y_t[i], s=300)
 plt.grid()

显示结果(绿色为离心力,红色为万有引力,蓝色为速度):

v_{y0} = 0.7·v_0:


v_{y0} = 2·v_0(超过第一宇宙速度,脱离该星体):

v_{y0} = v_0(做匀速圆周运动):

射电天文光谱数据处理3

本文中用到的实验数据是以 .fits 格式储存的。关于 .fits 文件格式,下边给出一些介绍:

FITS 全称是Flexible Image Transport System,它是天文学界常用的数据格式,它专门为在不同平台之间交换数据而设计,国际天文学会(IAU)于1982 年确定为世界各天文台之间用于数据传输、交换的统一标准格式。

例如实验用例:

scu = sc.read('orikl_hcop.fits',hdu=0)
print(scu)

其结果为

SpectralCube with shape=(7416, 105, 38) and unit=K:
 n_x:     38  type_x: RA---TAN  unit_x: deg    range:    83.746662 deg:   83.889447 deg
 n_y:    105  type_y: DEC--TAN  unit_y: deg    range:    -5.459508 deg:   -5.236409 deg
 n_s:   7416  type_s: FREQ      unit_s: GHz    range:      356.588 GHz:     356.814 GHz

 7416:  在频率(波长)维度有7416个通道
 105:   y方向像素点个数
 38:    x方向像素点个数
 unit_: 单位

展示出了一个天文模型在相对位置的辐射强度。

由于可能出现边缘数据缺省或数据错误,因此进行数据清洗、移除坏点的操作:

sc_data2 = np.nan_to_num(sc_data) # 移除坏点(变边缘处可能出现数据错误的情况,均赋值为0)
scu2 = sc(data=sc_data2,wcs=w0)

为了进行之后的可视化,在这里将频率轴的单位由 GHz 转换成 km/s,更能表现物理性质:

scu2b = scu2.with_spectral_unit(u.km / u.s, rest_value=356.73424 * u.GHz,  velocity_convention='radio')  
# 多普勒变换                                                                                                    
# 其中 rest_value 设置了一个于静止的频率                                                               
# 其余变换要相对于该静止频率变换
# 一般是一条谱线的跃迁值

展示此时根据数据获得的图像:


可以看到,此时图像里有两处非常明显的辐射。但是周围的地方由于辐射强度过弱并且区分不明显,因此显示不清楚。

接下来的操作是将所需要的图像部分进行对数处理,方便观察更细致的辐射图像。

im = ax.imshow(scm0[10:80,:], origin='lower', cmap='jet',norm=colors.LogNorm(vmin=vm0, vmax=vm1))  
# 将取vm0~vm1频率范围内的辐射进行对数处理

展示此时根据数据获得的图像(白色线条为经纬线,横坐标为赤经,纵坐标为赤纬):


学习心得

Python 语法较C,C++ 简单易懂,更适合用于数据处理、科学计算、Web 开发等领域。其可移植性和可扩展性都很强,并且有很强大的类库,能满足各种要求。在数据处理可视化方面,Python 更是拥有很大优势。

学习 Python 是一段有趣而充实的旅程。通过不断学习和实践,我建立了坚实的编程基础,并且对 Python 的应用领域有了更深入的了解。我相信,只要保持热情和持续学习,我会在 Python 编程领域不断成长。

参考文献

[1] 数据来源:国家天文科学数据中心

[2] 张韵华 王新茂 陈效群 张瑞 编 数值计算方法与算法(第三版) 科学出版社

[3] 天文数据基础与Python天文技术培训(2020-12) (china-vo.org)

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 213,186评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,858评论 3 387
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,620评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,888评论 1 285
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,009评论 6 385
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,149评论 1 291
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,204评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,956评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,385评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,698评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,863评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,544评论 4 335
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,185评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,899评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,141评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,684评论 2 362
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,750评论 2 351

推荐阅读更多精彩内容