今晚我在网上冲浪的时候看到了有老外写了怎么使用使用 Python 和 OpenCV 处理天文图像。大致浏览了一下，发现其中的一些原理其实也能用在遥感图像处理上。下面是我机翻的内容，分享出来给大家一起观摩、学习。以后你也可以变身孙连城，在仰望天空的时候，别忘了小白分享过这篇文章。--------------------------我是分割线--------------------------可能很多人不知道，天文相机的拍摄方式与普通相机截然不同。
我们不能随便拿起相机就去拍摄遥远的星系。遥远的深空天体，例如其他星系或星云，都非常暗淡，相机的矩阵灵敏度不足以捕捉到它们。
然而，传感器的噪声是随机的，如果我们把多张图像叠加起来，就可以提高信噪比，从而获得更高质量的图像。
在本教程中，我将使用放置在我阳台上的智能望远镜拍摄的图像。
我用这台望远镜拍摄了考德威尔19号星云（IC 5146）的图像，总曝光时间为一小时。我将使用Python和OpenCV处理这些图像，并展示其工作原理。
我们开始吧！
1. 硬件
一般来说，天文摄影是摄影中最具挑战性的部分。深空天体非常暗淡，为了收集足够的光线，通常需要大口径望远镜。而且我们都知道，地球也在自转，因此望远镜需要特殊的赤道仪来补偿这种自转。实际上，一套能够进行天文摄影的望远镜通常是一套重达20-30公斤的设备，例如：
我们还需要一台相机、一台笔记本电脑、一台用于自动导星相机的第二台望远镜、所有这些设备的外部电源等等。
听起来很复杂，实际上也确实如此。然而，近几年出现了一种新型望远镜，即所谓的智能望远镜。实际上，它并非“真正的”望远镜；它没有目镜，没有屏幕，也不能直接通过它进行观测。智能望远镜更像是便携式天文相机。更小的矩阵使得制造商可以使用相对较小的复消色差透镜，而缺少赤道仪则可以通过软件处理来弥补。因此，我们得到的是一个体积小得多的设备，虽然比普通相机略大，但足够紧凑（重量约为2.5公斤），可以放在阳台上或装进包里随身携带。
2. 加载数据
从最后一张照片可以看出，这台望远镜没有赤道仪。赤道仪可以让望远镜跟随地球自转：
这造成了一个主要限制——望远镜可以拍摄一系列图像，但每张图像的曝光时间限制为 10 秒，否则星星的轨迹就会显现出来。
如前所述，“智能望远镜”没有屏幕，它通过手机应用程序控制。该应用程序会自动进行图像处理，但我们也可以在应用程序设置中启用保存原始数据的功能。之后，我们可以下载所有图像，并获得 1-2 GB 的文件，如下所示：
我们可以看到，望远镜以所谓的FITS格式保存帧图像。这是一种在天文摄影中常用的格式，我们可以使用开源的Astropy库在Python中加载数据。import globfiles = glob.glob("IC 5146_sub/*.fit")files.sort()print("Files:", len(files))# > Files: 391print(files[:2])# > ['IC 5146_sub/Light_IC 5146_10.0s_LP_20251219-232930.fit',# > 'IC 5146_sub/Light_IC 5146_10.0s_LP_20251219-232941.fit']
借助 Astropy 库，我只需三行 Python 代码即可加载 FITS 图像。我们来试试第一张图像：
from astropy.io import fitsdef load_fits_data(filename: str) -> np.array: """ Load data from the FITS file """ with fits.open(filename) as image: return image[0].dataimage_data = load_fits_data(files[0])print("Type:", type(image_data))#> Type: <class 'numpy.ndarray'>print("Shape:", image_data.shape)#> Shape: (1920, 1080)
这里，我们有一个 NumPy 数组，1920x1080 是望远镜矩阵的分辨率。我可以使用 Matplotlib 绘制图像：
from matplotlib import pyplot as pltimport cv2plt.figure(figsize = (10, 4))image_rotated = cv2.rotate(image_data, cv2.ROTATE_90_CLOCKWISE)plt.imshow(image_rotated, cmap="gray", norm=LogNorm())plt.show()
结果如下，我们可以看到第一批星星：
原图是竖版的，我还旋转了一下以适应屏幕。这张照片目前还不太吸引人，肯定需要改进。让我们开始吧！
3. 预处理
首先，正如我们之前看到的，我们有一个 (1920, 1080) 形状的单色图像阵列，我们需要将其转换为 RGB 格式。大多数数码相机的矩阵使用微滤波器阵列，即所谓的拜耳滤波器，由布鲁斯·拜耳于 1976 年获得专利：
实际上，我们有一个包含不同颜色值的数组，我们可以将这些数据通过数学方法转换为“标准”的 RGB 格式。这个算法可能比较复杂，但我们可以使用 OpenCV 来实现。不过在此之前，让我们先绘制图像的直方图，看看亮度分布情况：
image_data = load_fits_data(files[0])print("Min:", np.min(image_data.flatten()))print("Max:", np.max(image_data.flatten()))#> Min: 0#> Max: 31662plt.figure(figsize =(10, 4))plt.gca().spines["top"].set_color('#dddddd')plt.gca().spines["right"].set_color('#dddddd')_ = plt.hist(image_data.flatten(), bins="auto", range=(0, 1000))
输出结果如下：
显然，原始数据是16位的（最大值为31662），但是大多数像素的值都在0到500之间。这很合理，因为这张天空图像是在夜晚拍摄的，当时天空很暗。
实际上，许多 OpenCV 算法都使用标准的 RGB 色彩空间，其中每个值都在 0 到 255 的范围内。考虑到 0 到 500 的范围，我可以将所有值除以 2 而不会丢失太多数据。现在，让我们将 RAW 图像转换为 RGB 图像：
k_reduce = 2def process_image(image: np.array) -> np.array: """ Convert a bayer matrix raw data to RGB """ data_out = np.uint8(np.clip(image / k_reduce, a_min=0, a_max=255)) rgb = cv2.cvtColor(data_out, cv2.COLOR_BAYER_GRBG2BGR) return rgb.astype(np.uint8)rgb = process_image(image_data)plt.figure(figsize = (10, 4))plt.imshow(cv2.rotate(rgb, cv2.ROTATE_90_CLOCKWISE))plt.show()
结果看起来好多了：
现在，我们可以看到一张可以处理的“正常”天空图像。细心的读者可能已经发现了中心的星云，虽然目前还不太清晰，但我们会努力改善它的效果。
4. 目标跟踪
正如开头提到的，图像堆叠的基本原理很简单。相机噪声是随机的，但物体是真实的。如果我们把多张图像组合起来，就能显著提高信噪比。一种“简单”的方法可以用几行代码实现，但我们先来看看图像效果。我随机选取了两张图像，并用 Matplotlib 将它们显示在一起：
img_rgb_1 = process_image(load_fits_data(files[0]))img_rgb_2 = process_image(load_fits_data(files[10]))f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 20))ax1.imshow(img_rgb_1)ax2.imshow(img_rgb_2)plt.show()
输出结果如下：
正如我们所见，“简单”的叠加方法行不通——所有恒星的位置都不一样。当然，如果所有图像都完美对齐，那么在一台价值3万5千元、配备 SkyWatcher/Takahashi 赤道仪和独立音频导星相机的望远镜上，这种方法或许可行。但对于入门级设备来说，情况并非如此。
为了对齐图像，我决定使用 OpenCV 的目标跟踪。为了提高跟踪精度，我会对图像进行预处理，提高亮度并将其转换为灰度图像。我还创建了一个方法来查找帧中最亮的星体，该星体将用作对齐的锚点：
def gamma(image: np.array, gamma: float) -> np.array: """ Apply gamma correction to the image """ lookup_table = np.empty((1,256), np.uint8) for i in range(256): lookup_table[0, i] = np.clip(pow(i / 255.0, gamma) * 255.0, 0, 255) return cv2.LUT(image, lookup_table)def preprocess(image: np.array) -> np.array: """ Preprocess image for the object tracking """ # Apply Gamma img = gamma(image, gamma=2.5) # Convert to grayscale img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) img = cv2.cvtColor(img, cv2.COLOR_GRAY2BGR) return np.clip(img, 0, 255).astype(np.uint8)def get_max_loc(image: np.array, bbox: List) -> Tuple[int, int]: """ Get the brightest pixel in the area """ (x, y, w, h) = [int(v) for v in bbox] crop_img = image[y:y+h, x:x+w] crop_gray = cv2.cvtColor(crop_img, cv2.COLOR_BGR2GRAY) crop_blurred = cv2.blur(crop_gray, (5,5)) minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(crop_blurred) return x + maxLoc[0], y + maxLoc[1]
现在，我可以手动选择两个具有可识别星形图案的区域：
bbox0 = (325, 319, 232, 232)bbox1 = (671, 1366, 238, 238)frame0 = process_image(load_fits_data(files[0]))frame0 = preprocess(frame0)def draw_on_frame(frame: np.array, bbox0: List, pt0: Tuple, bbox1: List, pt1: Tuple): """ Draw debug data on frame """ # Frame-1 (x, y, w, h) = [int(v) for v in bbox0] cv2.rectangle(frame, (x, y), (x+w, y+h), (200, 0, 0), 3) cv2.circle(frame, pt0, 7, (200, 0, 0), 3) # Frame-2 (x, y, w, h) = [int(v) for v in bbox1] cv2.rectangle(frame, (x, y), (x+w, y+h), (0, 200, 0), 3) cv2.circle(frame, pt1, 7, (0, 200, 0), 3)draw_on_frame(frame0, bbox0, pt0, bbox1, pt1)plt.figure(figsize = (16, 9))plt.imshow(frame0)plt.show()
用于调试的方法draw_on_frame是为了检查输出是否正确。结果如下所示：
思路很简单。我以第一张图像和两个感兴趣区域（ROI）作为参考。然后，对于其他每张图像，我使用目标跟踪找到相同的区域。最后，我可以将这些图像与第一张图像对齐。
在测试目标跟踪算法后，发现跟踪本身并非像素级精确。因此，我使用该区域内最亮的恒星作为锚点。现在让我们对所有图像运行目标跟踪：
frames = []frames_orig = []anchor_points = []# ROI (manually selected)frame0 = preprocess(images_rgb[0])bbox0 = (325, 319, 232, 232)bbox1 = (671, 1366, 238, 238)# Process all framesfor frame in tqdm(images_rgb[1:], desc="Processing"): tracker1 = cv2.legacy.TrackerBoosting_create() tracker1.init(frame0, bbox0) tracker2 = cv2.legacy.TrackerBoosting_create() tracker2.init(frame0, bbox1) frame_processed = preprocess(frame) ok0, bbox_0 = tracker1.update(frame_processed) ok1, bbox_1 = tracker2.update(frame_processed) if ok0 and ok1: pt0 = get_max_loc(frame, bbox_0) pt1 = get_max_loc(frame, bbox_1) anchor_points.append((pt0, pt1)) frames_orig.append(frame) frames.append(frame_processed)
至此，我们已经获得了所有帧的ROI坐标。实际上，它看起来是这样的：
现在我们可以进行配准了。
5. 图像对齐
有了数据之后，我们来校准一下所有星星！我会以第一帧为参考，移动并旋转其他帧。为什么要旋转呢？原来星星的位置不仅不同，而且由于地球自转，图像本身也会旋转——我们可以通过比较最后一张照片中红色和绿色矩形之间的角度来观察这一点。
当我们掌握了所有数据后，这个过程就很简单了：
# Frame-0f0_pt0, f0_pt1 = anchor_points[0]a0 = math.atan2(f0_pt1[1] - f0_pt0[1], f0_pt1[0] - f0_pt0[0])frames_processed = [ frames_orig[0] ]# Rotate and shift other frames to the first onefor ind in tqdm(range(1, len(frames)), desc="Processing"): frame = frames_orig[ind].copy() pt0, pt1 = anchor_points[ind] angle = math.atan2(pt1[1] - pt0[1], pt1[0] - pt0[0]) angle_diff_deg = math.degrees(angle - a0) # Rotate center = pt0 rotation_matrix = cv2.getRotationMatrix2D(center, angle_diff_deg, scale=1) frame = cv2.warpAffine(frame, rotation_matrix, (num_cols, num_rows)) # Shift dx1, dy1 = pt0[0] - f0_pt0[0], pt0[1] - f0_pt0[1] translation_matrix = np.float32([ [1,0,-dx1], [0,1,-dy1] ]) frame = cv2.warpAffine(frame, translation_matrix, (num_cols, num_rows)) # Save frames_processed.append(frame)
我们可以看到，图像被旋转和平移了，所有的星星现在都位于同一位置。6. 堆叠
现在，我们准备见证信号处理的神奇之处 :) 让我们把所有图像合并起来：
final_b = np.zeros(shape=(shape[0], shape[1]), dtype=np.float32)final_g = np.zeros(shape=(shape[0], shape[1]), dtype=np.float32)final_r = np.zeros(shape=(shape[0], shape[1]), dtype=np.float32)for image in tqdm(frames_rgb, desc="Processing"): # Split image to channels r, g, b = image[:, :, 0], image[:, :, 1], image[:, :, 2] final_r += r final_g += g final_b += bfinal_b = (final_b / count).astype(np.uint8)final_g = (final_g / count).astype(np.uint8)final_r = (final_r / count).astype(np.uint8)
由于我们使用的是原始图像，因此还可以应用白平衡。这里，我手动提高了红色通道的亮度：
k_red, k_green, k_blue = 1.5, 1.0, 1.0final_r = gamma(final_r, gamma=k_red)final_g = gamma(final_g, gamma=k_green)final_b = gamma(final_b, gamma=k_blue)
将各图层堆叠起来后，我们就可以创建最终图像。我还根据需要调整了亮度和对比度：
# Merge channels (BGR)final_image = cv2.merge([final_b, final_g, final_r])# Adjust contrast (alpha) and brightness (beta)final_image = cv2.convertScaleAbs(final_image, alpha=3.5, beta=-130)
效果令人惊艳。单帧画面如下所示：
最终，将 310 帧图像合成后，结果如下所示：
显然，在 Photoshop 或其他编辑软件中可以进行更多调整（色阶、曲线、锐化），但这超出了本文的范围。
感谢阅读。
10个月宝宝每天需要喝多少奶粉？