PixInsight、RC Astro 三件套、Siril 与 DeepSkyStacker：流程、数学模型与参考伪代码

目录

0. 边界、公开程度与阅读方式

1. 一张深空图像到底在解什么问题

深空后期不是“把照片修漂亮”这么简单。它本质上是在估计一个被光学系统、传感器、天空背景、热噪声、读出噪声、大气视宁度、导星误差、采样不足、光害梯度和离群污染共同破坏的天空亮度场。

其中 $R_k$ 是第 $k$ 张 light frame，$S$ 是天球坐标与波长上的真实信号。实际成片还要经过审美映射：人眼看不见线性数据中的微弱星云，所以必须做非线性拉伸、局部对比、颜色映射与星点控制。

2. 总流程图：从 RAW/FITS 到成片

RAW/FITS lights + bias/darks/flats/dark-flats
        │
        ├─► 1. Master calibration frames
        │       master_bias, master_dark, master_flat
        │
        ├─► 2. Calibrate each light
        │       subtract bias/dark, divide normalized flat, cosmetic correction
        │
        ├─► 3. Debayer or CFA drizzle preparation (OSC only)
        │
        ├─► 4. Star detection + registration
        │       select reference, estimate transforms, resample or drizzle metadata
        │
        ├─► 5. Normalize + weight + integrate
        │       local/global normalization, pixel rejection, weighted average
        │
        ├─► 6. Linear master
        │       crop, background model, color calibration, BXT, NXT
        │
        ├─► 7. Stretch
        │       MTF / Histogram / GHS / arcsinh
        │
        ├─► 8. Star separation / starless processing
        │       SXT, curves, contrast, saturation, star reduction
        │
        └─► 9. Final recombination and export
                TIFF/JPEG/PNG, optional Photoshop/GIMP finishing

PixInsight 的优势是几乎每个环节都有高级参数、可重复的 process icon 和脚本化流程；Siril 的优势是开源、脚本化、速度快、核心算法透明；DSS 的优势是低门槛和“把一堆 light/dark/flat 变成一个可继续后期的中间文件”。

3. 传感器与光学成像数学模型

3.1 光学卷积模型

若先忽略采样和传感器响应，光学系统把真实场景 $f$ 与点扩散函数 PSF $g$ 做卷积：

有噪声时：

这正是反卷积问题的起点。RC Astro 对 BXT 的数学说明也采用这个形式，并强调真实数据中 PSF 不准确、噪声存在、数值精度有限，所以“完美反卷积”不可达。

3.2 传感器 ADU 模型

对第 $k$ 张子帧、像素 $p$，可用下面的简化模型：

含义：

3.3 噪声方差

以电子数为单位，单帧像素的近似方差为：

如果转换到 ADU，方差会被增益平方缩放。这个模型解释了为什么深空摄影里“多拍、抖动、叠加”比后期插件更基础：后期可以改变噪声外观，但不能凭空增加信号中的信息量。

3.4 为什么必须保持线性

校准、反卷积、光度颜色校准等都默认像素值与入射光子数近似成比例。若先做非线性拉伸 $T$，则：

所以校准、配准、叠加、PCC/SPCC、BXT 这类步骤应在线性阶段完成。RC Astro 也明确要求 BXT AI4 使用线性输入，并不建议在 BXT 前做降噪。

4. 校准帧：Bias / Dark / Flat / Dark-flat

4.1 标准校准公式

最常见的 calibrated light 可写为：

其中 $M_D$ 是与 light 匹配温度、增益、offset、曝光时间的 master dark；$M_F$ 是 master flat；$M_{DF}$ 是 dark-flat 或 flat-dark；$\widetilde F$ 是归一化平场。

如果使用 bias 而不是 dark-flat，或 dark 需要按曝光时间缩放，可以写成更一般形式：

$\alpha_k=t_k/t_D$ 是暗场缩放因子。现代 CMOS 相机常不推荐任意缩放 dark，因为暗电流、amp glow 与传感器校正可能不严格线性；更稳妥的是拍匹配曝光的 dark 与 dark-flat。

4.2 参考伪代码：校准单张 light

function calibrate_light(light, master_dark, master_flat, master_darkflat=None):
    # All arrays should be linear, same size, same CFA pattern if OSC raw.
    if master_darkflat is not None:
        flat_signal = master_flat - master_darkflat
    else:
        flat_signal = master_flat - estimate_bias_or_offset(master_flat)

    flat_norm = flat_signal / robust_median(flat_signal)

    # Avoid division by zero or very weak flat pixels.
    flat_norm = clamp(flat_norm, lower=epsilon)

    calibrated = (light - master_dark) / flat_norm
    calibrated = cosmetic_correct_hot_cold_pixels(calibrated)
    return calibrated

4.3 校准帧的“减”和“除”分别在修什么

5. Master 帧构造、SNR 与误差传播

5.1 均值叠加的 SNR 推导

假设 $N$ 张独立子帧，每张同一像素的测量为：

均值为 $\bar X=\frac{1}{N}\sum_i X_i$，则：

这就是 Siril 文档里 average/sum stacking SNR 随 $\sqrt{N}$ 增长的数学来源。

5.2 中位数叠加的 SNR 近似

若噪声近似正态分布，中位数的渐近方差为：

正态分布在均值处的密度 $f(S)=1/(\sqrt{2\pi}\sigma)$，所以：

这也解释了为什么中位数常用于 master dark/flat/bias，而 light 更常用“均值 + 拒绝算法”：均值的 SNR 更高，拒绝算法负责剔除离群值。

5.3 加权均值的最优权重

若每张图像噪声方差不同，输出为：

在无偏约束下最小化方差，得到最优权重：

现实中除了背景噪声，还要考虑 FWHM、星点圆度、透明度、云、导星质量。Siril 提供按星数、weighted FWHM、背景噪声、积分时间等加权；PixInsight 1.8.9 引入了新的基于测光的图像质量估计和权重算法。

5.4 平场误差传播

校准公式 $C=(L-D)/F$ 中，$F$ 自身也有噪声。用一阶误差传播：

结论：平场不是越少越好；低 SNR 平场会把乘性噪声注入 light。平场 ADU 不应过低，数量也应足够。

5.5 参考伪代码：生成 master 帧

function make_master(frames, method="winsorized_sigma", normalize=None):
    # frames: list of same-shape linear arrays
    aligned = frames  # calibration frames usually do not require star alignment

    if normalize == "multiplicative":
        med_ref = median([robust_median(f) for f in aligned])
        aligned = [f * (med_ref / robust_median(f)) for f in aligned]

    master = empty_like(aligned[0])
    for pixel p:
        stack = [f[p] for f in aligned]
        kept = reject_outliers(stack, method=method)
        master[p] = mean(kept)  # or median for small/contaminated calibration sets
    return master

6. Debayer、CFA 与 Bayer Drizzle

6.1 Debayer 是插值，不是“恢复真实颜色”

OSC 相机每个像素只测一个颜色通道，例如 RGGB：

R G R G ...
G B G B ...
R G R G ...

普通 debayer 要从邻近像素推断缺失通道。双线性插值的参考形式：

更高级算法如 VNG/RCD 会利用边缘方向与空间/谱相关性减少彩边和伪色。Siril 文档也指出 debayer 算法依赖空间和光谱相关性去推断缺失分辨率。

6.2 Debayer 的伪代码

function bilinear_debayer(raw, pattern="RGGB"):
    R = zeros_like(raw); G = zeros_like(raw); B = zeros_like(raw)
    assign_known_CFA_samples(raw, R, G, B, pattern)

    for pixel p:
        if R[p] is missing: R[p] = average(nearest_known_R_neighbors(p))
        if G[p] is missing: G[p] = average(nearest_known_G_neighbors(p))
        if B[p] is missing: B[p] = average(nearest_known_B_neighbors(p))

    return RGB(R, G, B)

6.3 Bayer Drizzle 的思想

Bayer Drizzle 不先做颜色插值，而是把每个 CFA 原始像素作为“只属于某一通道的真实采样点”投到输出网格上。足够抖动和足够子帧数量时，输出的每个颜色通道都能得到更均匀采样。Siril 文档建议 OSC 常从 scale=1.0、pixfrac=1.0 开始；若 scale > 1.0，通道采样更稀疏，需要更多数据。

function cfa_drizzle(raw_frames, transforms, cfa_pattern, scale, pixfrac):
    out_R, weight_R = zeros(), zeros()
    out_G, weight_G = zeros(), zeros()
    out_B, weight_B = zeros(), zeros()

    for each frame k:
        T = transforms[k]  # raw pixel -> output coordinate
        for each raw pixel p:
            channel = cfa_color_at(p, cfa_pattern)  # R/G/B
            value = raw_frames[k][p]
            footprint = transformed_shrunken_pixel(p, T, scale, pixfrac)
            for output pixel q overlapping footprint:
                a = overlap_area(footprint, q)
                add value to corresponding channel:
                    out_C[q] += a * value
                    weight_C[q] += a

    return RGB(out_R/weight_R, out_G/weight_G, out_B/weight_B)

7. 星点检测、PSF 拟合与子帧质量指标

7.1 背景与噪声估计

星点检测的第一步通常是估计背景 $b$ 与噪声 $\sigma$。稳健估计常用：

$1.4826$ 是把 MAD 转成正态标准差的比例因子。

7.2 候选星点检测

Siril 的 Dynamic PSF 文档公开描述了流程：估计背景/噪声、用高斯核平滑、在背景加若干倍噪声阈值之上找局部最大、做 sanity check、识别饱和核心、再拟合 PSF 模型。

function detect_star_candidates(image, threshold_k, radius):
    b = median(image)
    sigma = 1.4826 * median(abs(image - b))
    smooth = gaussian_blur(image, sigma_guess)

    candidates = []
    for pixel p in smooth:
        if smooth[p] > b + threshold_k * sigma and is_local_maximum(smooth, p, radius):
            if core_and_neighbors_are_star_like(smooth, p, b, sigma):
                candidates.append(p)
    return candidates

7.3 PSF 模型：Gaussian 与 Moffat

星点可用二维高斯近似：

若 $\Sigma$ 的主轴标准差为 $\sigma_x,\sigma_y$，则：

真实星点常有更重的翼，可用 Moffat：

7.4 参考伪代码：PSF 拟合与质量指标

function fit_stars(image, candidates):
    stars = []
    for p in candidates:
        patch = crop(image, center=p, size=fit_box)
        params = nonlinear_least_squares(
            model="gaussian_or_moffat",
            data=patch,
            initial=estimate_from_derivatives(patch)
        )
        if params.amplitude > min_amp and params.fwhm in valid_range and params.eccentricity < max_e:
            stars.append(params)
    return stars

function frame_quality(stars, background_sigma):
    fwhm = median([s.fwhm for s in stars])
    ecc  = median([s.eccentricity for s in stars])
    nstars = len(stars)
    # Generic quality score, not a proprietary PI/Siril formula.
    score = nstars / (fwhm^2 * (1 + 4*ecc)^2 * background_sigma^2)
    return score

8. 配准：三角匹配、RANSAC、仿射/单应/畸变模型

8.1 为什么要配准

赤道仪误差、抖动、子午线翻转、dither、光学畸变和多晚拍摄会使同一颗星落在不同像素位置。叠加之前要估计每张图到参考图的几何变换 $T_k$。

8.2 常见变换模型

Siril 文档说明其全局配准基于三角形相似性匹配来识别共同星点，随后使用 RANSAC 排除 outlier 并确定投影矩阵；1.3 以后还可处理 SIP convention 的畸变，plate solving 可使用最高 5 阶多项式畸变。

8.3 仿射最小二乘推导

给定匹配星点 $(x_i,y_i)\to(u_i,v_i)$，仿射模型为：

写成矩阵：

其中 $W$ 可用星点测量误差或亮度构造。RANSAC 负责先找到可靠 inlier，最小二乘再在 inlier 上精修参数。

8.4 三角匹配 + RANSAC 伪代码

function register_frame(frame_stars, ref_stars, model="affine"):
    # 1. Build invariant descriptors from star triangles.
    tri_frame = build_triangles(frame_stars, max_triangles)
    tri_ref   = build_triangles(ref_stars, max_triangles)
    matches = match_triangles_by_side_ratios_and_angles(tri_frame, tri_ref)

    # 2. Convert triangle matches to candidate point correspondences.
    point_pairs = accumulate_star_pair_votes(matches)

    # 3. Robust model estimation.
    best_model, best_inliers = None, []
    for iter in 1..max_trials:
        sample = random_minimal_sample(point_pairs, model)
        T = fit_transform(sample, model)
        inliers = []
        for (x, y) in point_pairs:
            if distance(T(x), y) < tolerance_pixels:
                inliers.append((x, y))
        if len(inliers) > len(best_inliers):
            best_model, best_inliers = T, inliers

    # 4. Refit using all inliers.
    T_refined = weighted_least_squares_transform(best_inliers, model)
    return T_refined, best_inliers

8.5 重采样插值

配准后若直接保存 registered frames，要对输入图像做逆映射采样：

插值核可用 nearest、bilinear、bicubic、Lanczos 等。插值会引入相关噪声和轻微锐度损失；Drizzle 则尽量避免把输入像素先插值成输出像素。

9. 子帧归一化、局部归一化与权重

9.1 为什么 rejection 前要归一化

如果一张图天空背景高、另一张图背景低，同一像素 stack 的分布会变宽。Sigma clipping 会把正常背景差异误判为 outlier。Siril 文档明确说明 mean stacking with rejection 前的 normalization 很重要，因为月光、城市光害、温度变化等会造成图像 level 差异。

9.2 全局归一化公式

稳健背景和尺度：

几种常见归一化：

9.3 局部归一化参考模型

光害梯度、薄云、月光变化是空间变化的，不能只用一个全局 offset/scale。可设：

$a_i(p),b_i(p)$ 用网格采样后拟合平滑曲面。稳健目标函数：

$\rho$ 可用 Huber loss，$\lambda$ 控制曲面的平滑程度。

9.4 权重的一般形式

理论上单纯降低噪声应取 $w_i\propto1/\sigma_i^2$；成像质量还需要惩罚 FWHM 和星点拖线：

这不是 PixInsight 或 Siril 的精确公式，而是一个可解释的参考权重模型：透明度 $T_i$ 越高、背景噪声越低、FWHM 越小、偏心率越低，权重越高。

9.5 归一化与权重伪代码

function normalize_and_weight(frames, reference):
    stats_ref = robust_stats(reference)
    normalized = []
    weights = []

    for frame in frames:
        stats = robust_stats(frame)
        frame2 = stats_ref.median + (stats_ref.sigma / stats.sigma) * (frame - stats.median)

        stars = detect_and_fit_stars(frame2)
        fwhm = median(star.fwhm for star in stars)
        ecc = median(star.eccentricity for star in stars)
        bg_sigma = robust_background_sigma(frame2)
        transparency = estimate_transparency_from_star_fluxes(stars, reference)

        weight = (transparency^2 / bg_sigma^2) * (fwhm_ref / fwhm)^2 * exp(-4 * ecc)
        normalized.append(frame2)
        weights.append(weight)

    return normalized, weights

10. 叠加与像素拒绝：Average、Median、Sigma、MAD、Winsorized、LFC、GESD

10.1 加权叠加统一公式

给定已归一化、已配准图像 $I_i$，每个像素的输出为：

$m_i(p)\in\{0,1\}$ 是 rejection mask，表示该像素是否保留。

10.2 Sigma clipping

对某个像素位置的 stack $\{x_i\}$，迭代计算均值和标准差，剔除离均值超过阈值的点：

function sigma_clip(stack, k_low, k_high, max_iter):
    kept = stack
    for t in 1..max_iter:
        mu = mean(kept)
        sigma = std(kept)
        new_kept = [x for x in kept if mu-k_low*sigma <= x <= mu+k_high*sigma]
        if len(new_kept) == len(kept): break
        kept = new_kept
    return kept

10.3 MAD clipping

MAD clipping 用中位数和 MAD 代替均值/标准差，对小样本和强 outlier 更稳健：

function mad_clip(stack, k_low, k_high, max_iter):
    kept = stack
    for t in 1..max_iter:
        m = median(kept)
        sigma = 1.4826 * median(abs(kept - m))
        kept2 = [x for x in kept if m-k_low*sigma <= x <= m+k_high*sigma]
        if len(kept2) == len(kept): break
        kept = kept2
    return kept

10.4 Winsorized sigma clipping

Winsorization 的思想不是立即删除极端值，而是先把极端值压到边界，得到更稳健的 $\mu,\sigma$ 估计，再做 clipping。Siril 文档称其与 Sigma Clipping 类似，但对 outlier 检测更稳健，并参考 Huber robust statistics。

function winsorized_sigma_clip(stack, k_low, k_high, max_iter):
    kept = stack
    for t in 1..max_iter:
        med = median(kept)
        sigma0 = 1.4826 * median(abs(kept - med))
        lo = med - k_low * sigma0
        hi = med + k_high * sigma0
        wins = [min(max(x, lo), hi) for x in kept]
        mu = mean(wins)
        sigma = std(wins)
        kept2 = [x for x in kept if mu-k_low*sigma <= x <= mu+k_high*sigma]
        if len(kept2) == len(kept): break
        kept = kept2
    return kept

10.5 Median sigma clipping

Median Sigma Clipping 与 Sigma Clipping 的差别常在于被拒绝点不参与均值，而可用中位数替换，以避免空洞或极端值影响。

function median_sigma_integrate(stack, k_low, k_high):
    kept = sigma_clip(stack, k_low, k_high, max_iter)
    med = median(stack)
    repaired = [x if x in kept else med for x in stack]
    return mean(repaired)

10.6 Linear Fit Clipping：公开描述与参考实现

Siril 文档说明 Linear Fit Clipping 是 Juan Conejero / PixInsight ImageIntegration 相关算法，拟合像素 stack 的最佳直线 $y=ax+b$ 并拒绝 outlier，适合大样本、不同方向和分布的天空梯度、卫星/飞机轨迹。完整专有实现没有在本文复现；下面给出一个“公开描述等价思想”的参考实现。

function reference_linear_fit_clip(stack, k_low, k_high):
    # stack: normalized values x_i at one output pixel from N frames
    # Build robust expected positions using sorted Gaussian quantiles.
    xs = sort(stack)
    N = len(xs)
    q = [normal_quantile((r + 0.5) / N) for r in 0..N-1]

    # Fit xs ≈ a*q + b. For normal uncontaminated data, sorted values are nearly linear vs quantiles.
    a, b = robust_linear_regression(q, xs)  # Huber or Theil-Sen
    residuals = [xs[r] - (a*q[r] + b) for r in 0..N-1]
    sigma = 1.4826 * median(abs(residuals - median(residuals)))

    kept_sorted = []
    for r in 0..N-1:
        if -k_low*sigma <= residuals[r] <= k_high*sigma:
            kept_sorted.append(xs[r])
    return kept_sorted

10.7 Generalized ESD

Generalized Extreme Studentized Deviate Test 适合较大样本中一个或多个 outlier 检测。概念上每轮移除最极端点并计算统计量：

若 $R_j$ 超过临界值，则对应点可视为 outlier。Siril 文档称 GESD 对超过 50 张的大数据集表现很好。

10.8 不同 rejection 的适用性

11. Drizzle / CFA Drizzle 的数学与伪代码

11.1 Drizzle 的核心思想

Drizzle 由 Fruchter 与 Hook 提出，用于从欠采样、dithered 数据做线性重建。STScI 的介绍强调：输入像素会映射到输出子采样网格，考虑平移、旋转和畸变；为了避免再次用大像素 footprint 卷积图像，允许把输入像素收缩成 “drop”。

11.2 更新公式

设某个输入 drop 的像素值为 $i_{xy}$，权重为 $w_{xy}$，它与输出像素的重叠面积比例为 $a_{xy}$。输出像素已有值 $I_{xy}$ 和权重 $W_{xy}$，则更新为：

这等价于按重叠面积和输入权重做在线加权平均。

11.3 pixfrac 的取舍

11.4 Drizzle 伪代码

function drizzle_integrate(frames, transforms, weights, scale=2.0, pixfrac=0.8):
    output = zeros(output_shape(scale))
    weight = zeros(output_shape(scale))

    for k in 1..N:
        for each input pixel p:
            val = frames[k][p]
            w   = weights[k] * inverse_variance_or_mask(k, p)
            if is_rejected_or_bad(k, p): continue

            # Transform input pixel center and footprint to output coordinates.
            footprint = map_pixel_footprint(p, transforms[k], scale)
            drop = shrink_footprint_about_center(footprint, pixfrac)

            for q in output_pixels_overlapping(drop):
                a = fractional_overlap_area(drop, q)
                output[q] = (output[q]*weight[q] + val*w*a) / (weight[q] + w*a)
                weight[q] += w*a

    return output, weight

11.5 Drizzle 与普通重采样叠加的差异

12. 背景建模、梯度去除与平场残差

12.1 背景模型

线性 master 可表示为：

$S(p)$ 是星云/星系/星点信号，$B(p)$ 是光害、月光、薄云、平场误差导致的慢变化背景。背景去除要估计 $B$，但不能把真实大尺度星云误当背景。

12.2 多项式 / RBF / 样条模型

或使用径向基函数：

DBE 类工具本质上是“采背景点 → 剔除星点/目标污染 → 拟合平滑背景面 → 减法或除法校正”。

12.3 参考伪代码：DBE/ABE 风格背景去除

function background_extraction(image, sample_grid, model="thin_plate_spline"):
    samples = []
    for box in sample_grid:
        patch = crop(image, box)
        if patch_contains_nebula_or_large_star(patch):
            continue
        value = sigma_clipped_median(patch)
        samples.append((box.center_x, box.center_y, value))

    # Robust fit avoids contamination from missed stars/nebulosity.
    B = robust_surface_fit(samples, model=model, smoothness=lambda)

    if gradient_is_additive:
        corrected = image - B + median(B)
    else:
        corrected = image / (B / median(B))
    return corrected, B

12.4 加性还是乘性

13. PCC / SPCC 色彩校准的数学框架

13.1 颜色校准也是拟合问题

在 RGB 图像中，对第 $s$ 颗参考星，测得通道通量：

目录/光谱合成给出期望相对颜色 $q_{s,c}$。寻找通道增益 $g_c$：

也可以在颜色比值空间拟合，例如 $R/G$、$B/G$：

13.2 PCC 与 SPCC 的差别

Siril 文档强调 PCC/SPCC 必须在线性图像上执行；Siril 的 SPCC 使用 Gaia DR3 的光谱数据，并可使用本地或远程 catalogue。

13.3 参考伪代码：SPCC/PCC 风格

function photometric_color_calibration(image_rgb, wcs, sensor_response, filters):
    stars_img = detect_and_fit_stars(image_rgb.luminance())
    stars_cat = query_catalog(wcs, field_of_view, limit_mag)
    matches = match_image_stars_to_catalog(stars_img, stars_cat, wcs)

    observations = []
    predictions = []
    for match in matches:
        flux_obs = aperture_photometry_rgb(image_rgb, match.image_position)
        flux_pred = synthetic_rgb_flux(match.catalog_spectrum_or_color,
                                       sensor_response, filters)
        if photometry_is_clean(flux_obs, match):
            observations.append(flux_obs)
            predictions.append(flux_pred)

    gains = robust_fit_diagonal_gains(observations, predictions,
                                      white_reference="average_spiral_galaxy_or_solar")
    calibrated = apply_channel_gains(image_rgb, gains)
    return calibrated, gains, fit_residuals

14. 线性到非线性：MTF、arcsinh、曲线与掩膜

14.1 为什么拉伸是不可逆审美映射

线性 master 中，大部分深空信号接近黑场。拉伸把暗部展开、亮部压缩，是非线性映射：

一旦拉伸，像素差值不再与光子通量成比例，许多物理/统计假设失效。

14.2 Midtones Transfer Function, MTF

一个常见 MTF 形式为：

$m<0.5$ 时暗部被提升，$m>0.5$ 时图像变暗。

14.3 arcsinh stretch

arcsinh stretch 常用于提升暗部同时保留较好的星色：

$a$ 越大，暗部拉伸越强。

14.4 掩膜的数学表达

很多局部处理都是掩膜混合：

$M(p)\in[0,1]$。星点掩膜、亮度掩膜、范围掩膜、星云掩膜都属于这个框架。

15. 反卷积：经典算法、BXT 公开数学与参考伪代码

15.1 反卷积的病态性

频域中，卷积变乘法：

朴素逆滤波：

当 $|H(\omega)|$ 很小，高频噪声被剧烈放大，所以反卷积必须正则化、阻尼或限制目标 PSF。

15.2 Wiener 反卷积

$K$ 近似噪声/信号功率比。$K$ 越大，恢复越保守，噪声放大越少。

15.3 Tikhonov 正则化

$L$ 可以是梯度或 Laplacian，$\lambda$ 控制锐化与噪声/振铃之间的平衡。

15.4 Richardson-Lucy 迭代

对泊松噪声，经典 Richardson-Lucy 更新为：

$h^\star$ 是 PSF 的翻转。RL 保持非负并常用于天文，但迭代过多会产生振铃、暗环和噪声结构。

function richardson_lucy(y, psf, iterations, mask=None, damping=None):
    x = max(y, epsilon)
    psf_flip = flip(psf)
    for t in 1..iterations:
        estimate = convolve(x, psf) + epsilon
        ratio = y / estimate
        correction = convolve(ratio, psf_flip)
        if damping is not None:
            correction = damp_small_scale_noise(correction, y, damping)
        if mask is not None:
            x = mask * x * correction + (1-mask) * x
        else:
            x = x * correction
        x = max(x, 0)
    return x

15.5 BXT 公开数学：目标不是“零宽 PSF”

RC Astro 的 BXT 数学说明把真实反卷积结果描述为“原始图像 $f$ 与一个更小的新 PSF $g'$ 卷积，再加噪声与误差”：

若机器学习算法为 $\mathcal{F}[x,\mathbf W]$，输入训练图像为 $h=f*g+n$，则可写为：

这表达了一个重要思想：BXT 的目标不是把所有星变成数学点，而是把当前 PSF 变成受控、更小、更少像差的新 PSF。

15.6 星点与非星体分量分开处理

RC Astro 公开说明进一步把场景分为星点与非星体：

这解释了 BXT 为什么可以给星点和非星体细节设置不同 deconvolution/sharpening 强度。它是审美上有用的自由度，但不应用作严格科学测光处理。

15.7 BXT 公开功能边界

15.8 BXT 风格参考伪代码：非公开实现的可复现替代

function bxt_style_reference_deconvolution(linear_rgb, params):
    assert is_linear(linear_rgb)

    # 1. Detect stars and fit local PSFs.
    lum = luminance(linear_rgb)
    candidates = detect_star_candidates(lum, threshold_k=5, radius=3)
    stars = fit_stars(lum, candidates)

    # 2. Build a spatially varying PSF field.
    # Each tile gets a PSF estimated from nearby stars; fallback to manual FWHM if sparse.
    psf_field = interpolate_local_psf(stars, tile_size=512,
                                      fallback_fwhm=params.manual_fwhm)

    # 3. Separate stellar/nonstellar masks.
    star_mask = build_star_mask(stars, halos=True, saturation=True)
    nonstellar_mask = 1 - soft_dilate(star_mask)

    # 4. Process overlapping tiles.
    output_tiles = []
    for tile in overlapping_tiles(linear_rgb, size=512, overlap=64):
        psf = psf_field.at(tile.center)

        # Optional: correct asymmetric PSF toward a rounder PSF first.
        corrected = regularized_deconv(tile.data, psf,
                                       target_psf=round_psf(psf, params.correct_only_radius),
                                       lambda=params.regularization)

        # Stellar and nonstellar target PSFs can differ.
        star_part = corrected * star_mask[tile]
        ns_part   = corrected * nonstellar_mask[tile]

        star_deconv = regularized_deconv(star_part, psf,
                                         target_psf=shrink_psf(psf, params.star_sharpen))
        ns_deconv   = regularized_deconv(ns_part, psf,
                                         target_psf=shrink_psf(psf, params.nonstellar_sharpen))

        tile_out = blend_by_masks(star_deconv, ns_deconv, star_mask[tile])
        output_tiles.append(tile_out)

    # 5. Blend tiles smoothly to avoid seams.
    out = overlap_add_with_cosine_windows(output_tiles)

    # 6. Preserve flux/headroom as much as possible.
    out = prevent_negative_values(out)
    out = optional_flux_rescale_per_star(out, linear_rgb, stars)
    return out

15.9 BXT 在流程中的推荐位置

Recommended PI + BXT order:
  calibrated/registered/integrated master
  -> crop if needed
  -> remove obvious gradients
  -> channel combination if mono
  -> optional BXT Correct Only before SPCC for aberration centering
  -> SPCC / color calibration
  -> BXT sharpening, still linear
  -> NXT or other denoise
  -> stretch
  -> SXT / starless processing

16. 降噪：噪声模型、NXT 公开参数与参考算法

16.1 噪声组成

深空图像噪声既有高频像素噪声，也有大尺度低频 blotch；既有亮度噪声，也有颜色噪声。可把图像写成：

其中 $n_{LF}$ 常来自校准残差、背景建模误差、少量子帧统计不足或非线性处理放大。

16.2 多尺度分解

降噪常在多尺度上进行：

$D_j$ 是第 $j$ 个尺度的细节层，$A_J$ 是大尺度残差。软阈值：

若 $T_j=k\sigma_j$，可对不同尺度采用不同强度。

16.3 Non-local means 参考公式

非局部均值用相似 patch 加权平均：

优点是能保留重复结构；缺点是慢，且天文图像中真实微弱结构与噪声很难区分。

16.4 NXT 公开参数含义

RC Astro 的 NXT 2/AI3 手册公开了几个核心控制：Denoise 表示去除噪声量，1.00 意味着试图移除全部噪声但可能过度平滑；Iterations 使用 successive approximation 逐步去噪，更多 iteration 有时可保留高噪区域细节但也可能产生伪影；它还允许强度/颜色噪声分离，以及高频/低频噪声分离，并用 HF/LF Scale 定义二者分界。

16.5 NXT 风格参考伪代码

function nxt_style_reference_denoise(image, denoise, iterations,
                                     color_sep=True, freq_sep=True, hf_lf_scale=6,
                                     denoise_color=0.8, denoise_lf=0.4):
    x = image
    for t in 1..iterations:
        step = denoise / iterations

        if color_sep:
            Y, C1, C2 = convert_to_luminance_chrominance(x)
        else:
            Y, C1, C2 = x, None, None

        if freq_sep:
            Y_lf = gaussian_blur(Y, sigma=hf_lf_scale)
            Y_hf = Y - Y_lf

            # D_theta can be a CNN/transformer denoiser trained on astro data;
            # for an open reference, replace it with wavelet/NLM/BM3D.
            Y_hf_clean = learned_or_wavelet_denoise(Y_hf, strength=step)
            Y_lf_clean = learned_or_wavelet_denoise(Y_lf, strength=denoise_lf*step)
            Y_clean = Y_hf_clean + Y_lf_clean
        else:
            Y_clean = learned_or_wavelet_denoise(Y, strength=step)

        if color_sep:
            C1_clean = learned_or_wavelet_denoise(C1, strength=denoise_color*step)
            C2_clean = learned_or_wavelet_denoise(C2, strength=denoise_color*step)
            clean = convert_from_luminance_chrominance(Y_clean, C1_clean, C2_clean)
        else:
            clean = Y_clean

        x = blend(x, clean, amount=step)
        x = preserve_large_scale_flux(x, image)
    return x

16.6 降噪时机

17. 去星/分星：SXT 公开说明与参考算法

17.1 星点分离的数学模型

天文图像可写成：

去星的目标不是简单删除亮点，而是估计：

在线性图像中，星层差值更有物理意义；非线性图像中，差值更偏向视觉合成层。

17.2 SXT 公开使用要点

RC Astro 的 SXT usage notes 说明：去星非常困难，不可能对所有图像 100% 完美；SXT 训练覆盖了从相机镜头到 JWST 的广泛仪器，但严重光学缺陷仍可能失败。PixInsight 版本建议尽早使用，理想是在 integration 后仍为线性数据时；线性图像会在内部自动执行简单 MTF stretch，再精确反转回线性状态。

17.3 传统去星参考算法

传统去星可用形态学和 inpainting：

function classical_star_removal(image):
    lum = luminance(image)
    small_scale = lum - morphological_opening(lum, structuring_element)
    star_mask = threshold(small_scale, k * robust_sigma(small_scale))
    star_mask = grow_mask_to_include_halos_and_spikes(star_mask, lum)

    starless = inpaint(image, mask=star_mask, method="patchmatch_or_poisson")
    stars = image - starless
    return starless, stars, star_mask

这种方法容易在亮星、星云细丝、星系核心附近产生洞、糊斑或误删结构。

17.4 SXT 风格神经网络参考伪代码

function sxt_style_reference_star_separation(image, linear=True):
    if linear:
        # SXT public notes state a simple MTF-like stretch can be applied internally and reversed.
        stretch_params = estimate_mtf_for_visibility(image)
        work = mtf_stretch(image, stretch_params)
    else:
        work = image

    # 1. Neural segmentation: probability of stellar pixels, halos, spikes.
    P_star = star_segmentation_network(work)
    mask = hysteresis_threshold(P_star, low=0.2, high=0.5)
    mask = refine_mask_with_brightness_and_psf(work, mask)

    # 2. Context-aware inpainting / starless reconstruction.
    starless_work = background_reconstruction_network(work, mask)
    starless_work = enforce_no_large_scale_color_shift(starless_work, work, mask)

    # 3. Reverse stretch if linear input.
    if linear:
        starless = inverse_mtf(starless_work, stretch_params)
    else:
        starless = starless_work

    # 4. Star layer in linear arithmetic.
    stars = image - starless
    stars = max(stars, 0)  # optional; avoid negative halos, but inspect flux conservation
    return starless, stars, mask

17.5 星点重组公式

线性重组：

非线性 screen 类重组：

Screen 更像视觉合成，容易保持亮星存在感；线性加法更适合保持强度关系，但要避免 clipping。

17.6 分星处理流程

Linear master after BXT/NXT
  -> stretch main image
  -> SXT on linear or early-stretched image depending workflow
  -> starless: curves, saturation, local contrast, denoise if needed
  -> stars: mild stretch, color protection, optional star reduction
  -> recombine: linear add or screen-like blend
  -> final crop/color/export

18. PixInsight、Siril、DeepSkyStacker 的算法定位对比

18.1 总体定位

18.2 算法维度对比

18.3 版本与生态状态

• PixInsight 官方 FAQ 将其描述为为 astrophotography 和 technical imaging 设计的高级图像处理平台。
• Siril 1.4.4 于 2026-06-17 发布，1.4.x 已包含强大的 astrometry、stacking、SPCC、Drizzle/CFA Drizzle 与脚本改进。
• RC Astro 2026-06-24 宣布 BXT/NXT/SXT stand-alone CLI beta，不再依赖 PixInsight 或 Photoshop，并已可集成到 Siril、Seti Astro Suite Pro。
• DSS GitHub 仓库说明 DSS 是 2006 年创建的 freeware，2018 年开源，采用 BSD 3-Clause；GitHub releases 显示 6.2.1 为 2026-06-17 的 latest release。

19. 实战流程：PI 三件套、Siril、DSS 三条路线

19.1 PixInsight + BXT/NXT/SXT 推荐主流程

1. WBPP
   - calibrate lights with matched dark/flat/dark-flat
   - debayer or keep CFA drizzle path
   - register
   - integrate with appropriate rejection
   - optional drizzle integration

2. Linear master cleanup
   - DynamicCrop / Crop
   - DBE/ABE or equivalent background extraction
   - channel combination if mono
   - optional BXT Correct Only if aberration correction before SPCC is beneficial
   - SPCC / color calibration

3. Linear restoration
   - BXT: moderate star and nonstellar sharpening
   - NXT: conservative denoise; avoid 100% plastic background

4. Stretch
   - STF preview -> HistogramTransformation / GHS / arcsinh

5. Star separation and nonlinear enhancement
   - SXT to starless + stars
   - starless: curves, contrast, saturation, masks
   - stars: protect color, mild star reduction if needed
   - recombine

6. Final
   - color balance, crop, annotation if desired, export

19.2 Siril 路线

1. Convert raw/ser/fits to Siril sequence
2. Preprocess with calibration frames
3. Register with global registration or astrometric registration
4. Stack using average with rejection + normalization + weights
5. Crop / background extraction
6. PCC/SPCC on linear image
7. Optional RC-Astro CLI scripts: BXT/NXT/SXT if installed/licensed
8. Stretch, color, final export

Siril 的优势是流程透明与脚本化。对于大量 smart telescope、OSC、宽场和需要自动化的场景，Siril 1.4.x 以后非常有竞争力。

19.3 DSS 路线

1. Load lights, darks, flats, bias/dark-flats
2. Register checked pictures
3. Choose score threshold / best percentage
4. Stack with Average or Kappa-Sigma / Median Kappa-Sigma
5. Save 32-bit TIFF/FITS without aggressive DSS stretch
6. Continue in PixInsight / Siril / Photoshop / GIMP

DSS 最好被看作“生成干净中间文件的预处理器”。后续的梯度去除、SPCC、反卷积、AI 降噪、去星和精细曲线通常需要其他软件。

20. 常见伪影、诊断公式与修正策略

20.1 一个实用强度准则

许多后期伪影来自“把估计问题推到数据 SNR 的极限外”。可用简单准则：

也就是说，BXT/NXT/SXT 都应该被看作“更好地估计已存在信息”的工具，而不是替代曝光时间、对焦、导星、校准和 dither 的工具。

21. 端到端参考伪代码

21.1 完整开放参考流程

function deep_sky_pipeline(lights, darks, flats, darkflats, camera_info, options):
    # ---------- Calibration masters ----------
    master_dark = make_master(darks, method="winsorized_sigma", normalize=None)
    master_darkflat = make_master(darkflats, method="winsorized_sigma", normalize=None)
    master_flat = make_master(flats, method="winsorized_sigma", normalize="multiplicative")

    # ---------- Calibrate lights ----------
    calibrated = []
    for L in lights:
        C = calibrate_light(L, master_dark, master_flat, master_darkflat)
        C = cosmetic_correct_hot_cold_pixels(C)
        calibrated.append(C)

    # ---------- Debayer or keep CFA path ----------
    if camera_info.is_OSC and not options.cfa_drizzle:
        calibrated = [debayer(C, method=options.debayer_method) for C in calibrated]

    # ---------- Star detection and reference selection ----------
    star_lists = [fit_stars(C, detect_star_candidates(luminance(C), 5, 3)) for C in calibrated]
    qualities = [frame_quality(stars, robust_background_sigma(luminance(C)))
                 for C, stars in zip(calibrated, star_lists)]
    ref_index = argmax(qualities)
    reference = calibrated[ref_index]

    # ---------- Registration ----------
    transforms = []
    for stars in star_lists:
        T, inliers = register_frame(stars, star_lists[ref_index], model=options.registration_model)
        transforms.append(T)

    # ---------- Integration ----------
    if options.cfa_drizzle:
        master_linear = cfa_drizzle(calibrated, transforms, camera_info.cfa_pattern,
                                    scale=options.drizzle_scale, pixfrac=options.pixfrac)
    elif options.drizzle:
        master_linear = drizzle_integrate(calibrated, transforms, weights=qualities,
                                          scale=options.drizzle_scale, pixfrac=options.pixfrac)
    else:
        registered = [resample(C, T, reference_grid=reference.grid) for C, T in zip(calibrated, transforms)]
        normalized, weights = normalize_and_weight(registered, reference)
        master_linear = integrate_with_rejection(normalized, weights,
                                                 rejection=options.rejection)

    # ---------- Linear processing ----------
    master_linear = crop_bad_edges(master_linear)
    master_linear, background_model = background_extraction(master_linear, auto_or_manual_samples())
    master_linear = photometric_color_calibration(master_linear, solve_wcs(master_linear),
                                                  camera_info.sensor, camera_info.filters)

    # BXT-style step: use real BXT if available, otherwise reference deconvolution.
    if options.use_bxt:
        master_linear = run_BXT(master_linear, options.bxt_params)
    else:
        master_linear = bxt_style_reference_deconvolution(master_linear, options.deconv_params)

    # NXT-style step.
    if options.use_nxt:
        master_linear = run_NXT(master_linear, options.nxt_params)
    else:
        master_linear = nxt_style_reference_denoise(master_linear, options.denoise)

    # ---------- Stretch and star workflow ----------
    stretched = stretch(master_linear, method=options.stretch_method)

    if options.use_sxt:
        starless, stars = run_SXT(stretched or master_linear, return_stars=True)
    else:
        starless, stars, mask = sxt_style_reference_star_separation(stretched, linear=False)

    starless = enhance_nebula_or_galaxy(starless, masks=build_masks(starless))
    stars = process_star_layer(stars, preserve_color=True, reduce_size=options.star_reduce)
    final = recombine_starless_and_stars(starless, stars, method=options.recombine)
    final = final_color_contrast_crop_export(final)
    return final

21.2 PixInsight 等价流程的“过程图标”思维

WBPP
  -> Blink / SubframeSelector
  -> ImageIntegration diagnostic review
  -> DynamicCrop
  -> DBE or ABE
  -> ChannelCombination if mono
  -> SPCC
  -> BlurXTerminator
  -> NoiseXTerminator
  -> HistogramTransformation or GHS
  -> StarXTerminator
  -> Curves / LocalHistogramEqualization / ColorSaturation on starless
  -> Curves / MorphologicalTransformation or star reduction on stars
  -> PixelMath recombination
  -> final Histogram / SCNR only if justified / export

21.3 Siril 脚本式思维

# Conceptual Siril-like pipeline, not exact command syntax for every version.
convert raw_files -> sequence
preprocess sequence -dark=master_dark -flat=master_flat -darkflat=master_darkflat
register sequence -method=global_or_astrometric -disto=on
stack sequence -method=average -rejection=winsorized_or_mad -norm=addscale -weight=wfwhm
crop
background_extraction
platesolve
spcc
# Optional from 2026 RC-Astro CLI integration:
run_rcastro_bxt
run_rcastro_nxt
stretch
run_rcastro_sxt or native star-removal script
final curves/export

21.4 DSS 最小高质量思路

load lights/darks/flats/bias_or_darkflats
check all
register pictures
inspect scores/FWHM-like quality, keep best 70-90% unless many bad frames
stack with:
    N < 10: average or median depending contamination
    N >= 10: kappa-sigma / median kappa-sigma
save 32-bit TIFF/FITS with embedded adjustments off or minimal
continue linear-ish processing in Siril/PixInsight when possible

22. 资料来源

以下来源用于核对版本、公开功能和算法说明。商业软件/插件未公开的内部实现，本文没有当作事实引用。

1. PixInsight FAQ：https://pixinsight.com/faq/
2. PixInsight ImageWeighting 文档：https://pixinsight.com/doc/docs/ImageWeighting/ImageWeighting.html
3. PixInsight SPCC 文档：https://pixinsight.com/doc/docs/SPCC/SPCC.html
4. RC Astro BlurXTerminator Technical Manual：https://www.rc-astro.com/blurxterminator-technical-manual/
5. RC Astro, The Mathematics of BlurXTerminator：https://www.rc-astro.com/the-mathematics-of-blurxterminator/
6. RC Astro NoiseXTerminator 2/AI3 User Manual：https://www.rc-astro.com/noisexterminator-2-ai3-user-manual-pixinsight/
7. RC Astro StarXTerminator Usage Notes：https://www.rc-astro.com/starxterminator-usage-notes/
8. RC Astro Stand-Alone Tools：https://www.rc-astro.com/stand-alone-rc-astro-tools/
9. Siril 1.4.4 Release：https://siril.org/download/2026-06-17-siril-1-4-4/
10. Siril Stacking documentation：https://siril.readthedocs.io/en/stable/preprocessing/stacking.html
11. Siril Registration documentation：https://siril.readthedocs.io/en/stable/preprocessing/registration.html
12. Siril Dynamic PSF documentation：https://siril.readthedocs.io/en/stable/Dynamic-PSF.html
13. Siril Drizzle documentation：https://siril.readthedocs.io/en/stable/preprocessing/drizzle.html
14. Siril SPCC documentation：https://siril.readthedocs.io/en/latest/processing/color-calibration/spcc.html
15. Siril Platesolving documentation：https://siril.readthedocs.io/en/stable/astrometry/platesolving.html
16. DeepSkyStacker GitHub repository：https://github.com/deepskystacker/DSS
17. DeepSkyStacker Technical Info：https://deepskystacker.free.fr/english/technical.htm
18. Fruchter & Hook, Drizzle: A Method for the Linear Reconstruction of Undersampled Images：https://arxiv.org/abs/astro-ph/9808087
19. STScI Drizzle introduction：https://www.stsci.edu/~fruchter/dither/drizzle.html
20. RANSAC documentation / Fischler & Bolles reference：https://www.nv5geospatialsoftware.com/docs/ransac.html
21. SExtractor paper / source extraction reference：https://ui.adsabs.harvard.edu/abs/1996A%26AS..117..393B/abstract