如何使用specutils/astropy将高分辨率1D光谱降级至目标分辨率？

核心原理：卷积的分辨率合成 #

光谱分辨率通常定义为 ( R = \lambda / \Delta\lambda )，其中 ( \Delta\lambda ) 是光谱特征的半高全宽（FWHM）。用高斯核平滑高分辨率光谱时，最终分辨率由模型固有分辨率和平滑核分辨率共同决定，遵循卷积的FWHM合成公式：
[
\text{exp_FWHM}^2 = \text{model_FWHM}^2 + \text{kernel_FWHM}^2
]
因此，所需平滑核的FWHM为：
[
\text{kernel_FWHM} = \sqrt{\text{exp_FWHM}^2 - \text{model_FWHM}^2}
]
若模型分辨率远高于实验分辨率（( \text{model}_R \gg \text{exp}_R )），可近似认为 ( \text{kernel_FWHM} \approx \text{exp_FWHM} )；但两者差距不大时，必须使用精确公式计算。

步骤1：将分辨率转换为FWHM #

已知模型分辨率 ( R_{\text{model}} ) 和实验分辨率 ( R_{\text{exp}} )：

• 对任意波长 ( \lambda )，模型固有FWHM：( \Delta\lambda_{\text{model}} = \lambda / R_{\text{model}} )
• 实验目标FWHM：( \Delta\lambda_{\text{exp}} = \lambda / R_{\text{exp}} )

如果是常数分辨率（R不随波长变化），可对全光谱用统一核；如果是随波长变化的分辨率，则需为每个波长点计算对应sigma。

步骤2：计算平滑核的sigma #

用你熟悉的公式将核的FWHM转换为高斯核标准差sigma：
[
\sigma = \text{kernel_FWHM} / (2 \times \sqrt{2 \times \ln2})
]
注意：sigma单位需与光谱波长轴单位一致（如埃、纳米），不能用像素单位！这是常见易错点。

步骤3：代码实现 #

情况1：使用specutils的gaussian_smooth #

import numpy as np
from specutils import Spectrum1D
from specutils.manipulation import gaussian_smooth

# 假设已加载模型光谱：model_spec（Spectrum1D对象）
model_R = 10000  # 模型分辨率
exp_R = 1000     # 实验分辨率

# 获取波长值（确保单位统一）
lambda_vals = model_spec.spectral_axis.value
# 计算各波长点的FWHM
model_FWHM = lambda_vals / model_R
exp_FWHM = lambda_vals / exp_R

# 计算核的FWHM（避免根号内为负，确保实验分辨率低于模型）
kernel_FWHM = np.sqrt(np.maximum(exp_FWHM**2 - model_FWHM**2, 0))
# 转换为sigma
sigma = kernel_FWHM / (2 * np.sqrt(2 * np.log(2)))

# 执行平滑（specutils支持随波长变化的sigma）
smoothed_spec = gaussian_smooth(model_spec, stddev=sigma)

# （可选）重采样到实验光谱的波长轴
from specutils.manipulation import FluxConservingResampler
from astropy import units as u

# 假设exp_wavelengths是实验的波长轴（带单位）
resampler = FluxConservingResampler()
resampled_spec = resampler(smoothed_spec, exp_wavelengths)

情况2：使用astropy的convolve #

astropy卷积默认用像素为单位，需将sigma转换为像素单位：

import numpy as np
from astropy.convolution import Gaussian1DKernel, convolve
from specutils import Spectrum1D

model_R = 10000
exp_R = 1000
model_spec = Spectrum1D(...)  # 你的模型光谱

lambda_vals = model_spec.spectral_axis.value
# 假设波长轴均匀采样，获取单像素对应的波长间隔
delta_lambda = np.diff(lambda_vals)[0]

# 计算核的FWHM（波长单位）
model_FWHM = lambda_vals / model_R
exp_FWHM = lambda_vals / exp_R
kernel_FWHM = np.sqrt(np.maximum(exp_FWHM**2 - model_FWHM**2, 0))

# 转换为像素单位的sigma
sigma_wavelength = kernel_FWHM / (2 * np.sqrt(2 * np.log(2)))
sigma_pixel = sigma_wavelength / delta_lambda

# 常数分辨率下取平均sigma，生成核
avg_sigma_pixel = np.mean(sigma_pixel)
kernel = Gaussian1DKernel(stddev=avg_sigma_pixel)

# 卷积通量并重构光谱对象
smoothed_flux = convolve(model_spec.flux.value, kernel)
smoothed_spec = Spectrum1D(flux=smoothed_flux * model_spec.flux.unit,
                           spectral_axis=model_spec.spectral_axis)

# （可选）重采样到实验波长轴
from specutils.manipulation import FluxConservingResampler
resampler = FluxConservingResampler()
resampled_spec = resampler(smoothed_spec, exp_wavelengths)

验证方法 #

确认平滑效果是否达标：

• 在模型光谱中插入一条已知FWHM的窄线（远小于模型固有FWHM）
• 平滑后测量该线的FWHM，看是否等于实验目标FWHM
• 若结果吻合，说明平滑参数设置正确

内容的提问来源于stack exchange，提问作者gabriel avellaneda