科学技术原理
Python 以及衍生的一系列数据处理分析软件包已经在天文领域获得了广泛使用,在天文数据处理和分析、数据管理、天文数据绘图、高性能计算、机器学习和深度学习、乃至望远镜管理等诸多领域,都有着广泛的使用和不俗的性能。这里着重介绍可视化的应用,即充分利用 matplotlib
库。本次任务会介绍相对简单的两种表现:第一,模拟天体运动轨道并可视化;第二,射电天文光谱数据可视化。
1. 天体运动可视化:通过物理公式计算模拟轨道的参数并将其显示在图像上
2. 射电天文光谱数据可视化:通过简单的数据提取和数据处理将结果显示在图像上1
设计方案
天体运动可视化
关于物理公式
A天体(质量为m)在与B天体(质量为M)作用下的万有引力的影响下按一定轨道运动。 如果我们用数学公式表示A天体此时接收到的力:
运动方程:
初始位置设置在x轴上,给一个y轴方向的初速度:
将运动方程改写成一阶微分方程(要使用使用四阶 Runge-Kutta
法求解):
为了在程序上便于操作,将变量归纳为矩阵形式
因此把一阶微分方程改写成如下形式:
关于算法(四阶Runge-Kutta法)2
四阶 Runge-Kutta
法用于求常微分方程的较高精度的数值解。算法表示如下:
其中O是零矩阵
初始值:
射电天文光谱数据处理
由于实验数据格式的特殊性,需要用到额外的库进行数据处理
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()
显示结果(绿色为离心力,红色为万有引力,蓝色为速度):
射电天文光谱数据处理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] 张韵华 王新茂 陈效群 张瑞 编 数值计算方法与算法(第三版) 科学出版社