Astropy实战：从FITS文件到科学级天文图像处理

1. 初识Astropy与FITS文件

第一次接触天文数据处理时，我被那些神秘的FITS文件搞得一头雾水。这种Flexible Image Transport System格式就像天文界的"黑匣子"，里面不仅存储着图像数据，还包含了观测时间、望远镜参数等关键元数据。记得当时我用普通图片查看器打开FITS文件，结果只看到一堆乱码，那种挫败感至今记忆犹新。

Astropy的出现彻底改变了这种局面。这个专门为天文学设计的Python库，就像是给天文数据研究者的一把瑞士军刀。通过几行简单的代码，我们就能轻松读取FITS文件：

from astropy.io import fits # 打开FITS文件 hdul = fits.open('observation.fits') # 查看文件结构 hdul.info()

运行这段代码，你会看到类似这样的输出：

Filename: observation.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 220 (2048, 2048) float32 1 SCI 1 ImageHDU 53 (2048, 2048) float32

这告诉我们这个FITS文件包含两个扩展：PRIMARY主头和SCI科学数据。每个扩展都有自己的一组头卡（Cards）和数据维度。对于刚入门的朋友，可以这样理解：FITS文件就像个多层蛋糕，每一层都存储着不同类型的信息。

2. 科学级数据清洗技巧

拿到原始FITS数据后，直接进行分析就像用未校准的望远镜观测——结果往往不可靠。我们需要先进行数据清洗，这包括去除坏像素、本底校正和平场处理等步骤。

2.1 处理坏像素和宇宙线

坏像素就像照片上的污点，会严重影响后续分析。Astropy提供了多种工具来识别和处理这些问题：

import numpy as np from astropy.stats import sigma_clip # 假设data是我们的科学数据 sci_data = hdul[1].data # 使用sigma clipping识别异常值 masked_data = sigma_clip(sci_data, sigma=3, maxiters=5) # 用中值替换异常值 median_val = np.ma.median(masked_data) cleaned_data = masked_data.filled(median_val)

这里sigma_clip函数会迭代地排除偏离中值3σ以上的像素，非常适合处理宇宙线击中产生的尖峰噪声。

2.2 本底和平场校正

本底噪声就像相片的暗电流，而平场校正则类似于消除镜头渐晕。Astropy配合Numpy可以轻松完成这些操作：

# 假设我们已经加载了bias和flat场文件 bias = fits.getdata('bias.fits') flat = fits.getdata('flat.fits') # 本底校正 corrected_data = (sci_data - bias) # 平场归一化 norm_flat = flat / np.median(flat) final_data = corrected_data / norm_flat

记得在处理平场时要先归一化，否则会改变数据的绝对亮度尺度。我曾经犯过这个错误，导致后续测光结果全部偏差了一个比例因子。

3. 天体坐标系的魔法转换

天文观测中最让人头疼的莫过于各种坐标系转换了。赤道坐标、银道坐标、地平坐标...就像不同的语言描述同一个天体。Astropy的coordinates模块让这些转换变得轻而易举。

3.1 基本坐标转换

假设我们从星表查到一个天体的赤道坐标，想转换为观测时的地平坐标：

from astropy.coordinates import SkyCoord, EarthLocation, AltAz from astropy.time import Time import astropy.units as u # 定义天体赤道坐标 icrs = SkyCoord(ra=146.3*u.deg, dec=-23.5*u.deg, frame='icrs') # 定义观测时间和地点 obs_time = Time('2023-06-15 02:00:00') keck = EarthLocation.of_site('Keck') # 自动获取凯克望远镜位置 # 转换为地平坐标系 altaz = icrs.transform_to(AltAz(obstime=obs_time, location=keck)) print(f"高度角: {altaz.alt:.2f}, 方位角: {altaz.az:.2f}")

这个功能在观测规划时特别有用。记得有一次我忘了考虑坐标系转换，结果望远镜指向完全错误的位置，白白浪费了宝贵的观测时间。

3.2 坐标匹配实战

当我们需要将不同星表中的天体进行匹配时，match_coordinates_sky函数就是救星：

from astropy.coordinates import match_coordinates_sky # 假设cat1和cat2是两个星表坐标 cat1 = SkyCoord(ra=[10.68458, 11.23125]*u.deg, dec=[41.26917, 42.12347]*u.deg) cat2 = SkyCoord(ra=[10.12345, 11.23120]*u.deg, dec=[42.24680, 42.12350]*u.deg) # 匹配两个星表 idx, d2d, _ = match_coordinates_sky(cat1, cat2) # 输出匹配结果 for i, (sep, matched_idx) in enumerate(zip(d2d, idx)): print(f"源{i}最佳匹配是{matched_idx}号源，角距离{sep.to(u.arcsec):.2f}")

在实际项目中，我通常设置一个合理的匹配半径（如1角秒），筛选出d2d小于该值的匹配对，这样可以有效避免误匹配。

4. 专业级天文图像可视化

原始的天文图像往往对比度很低，直接显示很难看出细节。Astropy的visualization模块提供了多种专业级的图像拉伸方法。

4.1 基本图像显示

让我们从最基本的显示开始：

import matplotlib.pyplot as plt from astropy.visualization import simple_norm # 使用ZScale算法自动调整显示范围 norm = simple_norm(final_data, 'zscale') plt.figure(figsize=(10, 8)) plt.imshow(final_data, norm=norm, origin='lower', cmap='gray') plt.colorbar(label='Flux (ADU)') plt.title('处理后的科学图像') plt.show()

这里的zscale算法是天文学家常用的自动对比度调整方法，它能很好地显示图像中的弱信号而不至于过曝亮源。

4.2 高级色彩映射

对于科学分析，选择合适的色彩映射至关重要。Astropy提供了多种专业拉伸方法：

from astropy.visualization import (ImageNormalize, LogStretch, AsinhStretch, ZScaleInterval) # 创建不同的标准化对象 norm1 = ImageNormalize(final_data, interval=ZScaleInterval(), stretch=LogStretch()) norm2 = ImageNormalize(final_data, interval=ZScaleInterval(), stretch=AsinhStretch(a=0.1)) # 绘制对比图 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6)) im1 = ax1.imshow(final_data, norm=norm1, origin='lower', cmap='viridis') ax1.set_title('对数拉伸') plt.colorbar(im1, ax=ax1) im2 = ax2.imshow(final_data, norm=norm2, origin='lower', cmap='plasma') ax2.set_title('反双曲正弦拉伸') plt.colorbar(im2, ax=ax2)

LogStretch适合显示动态范围极大的图像（如同时包含亮星和暗弱星系），而AsinhStretch则在保留亮源细节的同时展现弱信号。根据我的经验，处理星系图像时a=0.1是个不错的起点，而球状星团可能需要调整到0.05左右。

5. 从数据到科学图像的全流程实战

让我们把这些技术整合到一个完整的处理流程中。假设我们有一组原始观测数据，目标是生成一张可用于发表的科学图像。

5.1 数据准备与校准

from astropy.io import fits import numpy as np # 加载所有帧 science_files = ['science1.fits', 'science2.fits', 'science3.fits'] bias = fits.getdata('master_bias.fits') flat = fits.getdata('master_flat.fits') # 处理每一帧 processed = [] for file in science_files: with fits.open(file) as hdul: data = hdul[0].data # 数据校准 calibrated = (data - bias) / flat # 去除宇宙线 masked = sigma_clip(calibrated, sigma=3, maxiters=5) cleaned = masked.filled(np.ma.median(masked)) processed.append(cleaned) # 中值叠加 final_image = np.median(processed, axis=0)

这个流程中，中值叠加（median combine）能有效消除瞬态噪声（如宇宙线），同时保留真实天体信号。对于深度观测，我通常会使用20帧以上的叠加。

5.2 WCS坐标系统处理

科学图像必须包含准确的天球坐标信息。Astropy能完美处理FITS头中的WCS信息：

from astropy.wcs import WCS # 从原始文件复制WCS信息 with fits.open(science_files[0]) as hdul: wcs = WCS(hdul[0].header) # 保存最终图像 hdu = fits.PrimaryHDU(final_image, header=wcs.to_header()) hdu.writeto('final_image.fits', overwrite=True)

5.3 科学级图像输出

最后，我们需要生成适合出版的图像：

from astropy.visualization import make_lupton_rgb # 假设我们有g、r、i三个波段的图像 r_band = fits.getdata('final_r.fits') g_band = fits.getdata('final_g.fits') i_band = fits.getdata('final_i.fits') # 创建RGB图像 rgb = make_lupton_rgb(i_band, r_band, g_band, stretch=0.5, Q=10) plt.figure(figsize=(12, 10)) plt.imshow(rgb, origin='lower') plt.axis('off') plt.savefig('science_rgb.png', bbox_inches='tight', dpi=300)

这里的make_lupton_rgb采用了Lupton(2004)提出的算法，能在不丢失科学信息的前提下生成美观的彩色图像。stretch参数控制对比度，Q参数调整非线性程度，通常需要根据图像特点反复调试。