RedshiftExtrapolator 使用与原理说明#

本文档详细解释 autohea.core.utils.RedshiftExtrapolator 的物理推导、计算流程与使用方法，重点阐述“在给定信噪比阈值下可探测的最大红移”的计算原理。

源码位置：autohea/src/autohea/core/utils.py
相关函数：compute()、find_redshift_for_snr()、snr_li_ma()

1. 问题背景与目标#

目标：在给定的XSPEC物理模型、仪器响应（ARF/RMF）、背景与曝光时间条件下，求在某个SNR阈值（如 7σ）时的最大可探测红移 z_max。

核心变量：

模型 M(E; params)：XSPEC光子数谱（photons cm⁻² keV⁻¹ s⁻¹）
仪器响应：ARF（有效面积）与 RMF（能量重分配）
背景计数与时长：bkgnum, duration
SNR阈值：snr_target（默认7）

2. 物理论证与关键公式#

记观测到的光子数谱为 N_obs(E_obs)。类中遵循以下“严格的XSPEC红移外推公式”：

几何衰减（真实物理距离）：使用共动距离 r_c(z)

光子数密度按真实距离平方反比衰减：∝ 1/r_c²
因此归一化的几何因子：

r_c 因子 = [r_c(z₀) / r_c(z)]²

K-correction（幂律谱示例）

若 N_rest(E) ∝ E^(-α)，则红移导致能量变换 E_rest = E_obs × (1+z)
XSPEC的单位与红移下时间膨胀、能量间隔效应相互抵消，最终只需要能谱幂律的额外修正：

K 因子 = ((1+z₀) / (1+z))^α

归一化缩放总因子

设 norm_original 为 z₀ 时模型的归一化，z 时应使用：

norm_new = norm_original × ((1+z₀)/(1+z))^α × [r_c²(z₀)/r_c²(z)]

这正是代码中 find_redshift_for_snr() 计算 total_factor 的来源：

geometric_factor = (r_c_z0 / r_c_grid) ** 2
k_correction_factor = ((1 + z0) / (1 + z_grid)) ** alpha
total_factor = k_correction_factor * geometric_factor

注：光度距离 D_L = (1+z) r_c 是定义量，适合用于能量通量 F 与光度 L 的关系；但对光子数密度与XSPEC模型的处理，采用真实距离 r_c 更能直接体现几何衰减，因此本类选择 r_c。

3. 计算流程（逐步说明）#

以 compute(snr_target=7) → find_redshift_for_snr() 的实际实现为主线：

3.1 模型初始化与参数注入#

init_model() 内部依次执行：
- 初始化HEASoft环境，创建 xspec.Model(self._model) 实例；
- _set_par() 将传入的 par 按模型参数顺序写入XSPEC，并冻结参数；
- 自动识别红移参数（若存在 ...Redshift）与归一化参数（通常为最后一组件的 norm）。

3.2 红移网格与缩放因子#

初始红移区间 [zmin, zmax]（默认 [z0, z0+1]），生成均匀 z_grid；
依次计算：
- r_c_z0 = comoving_distance(z0)；
- r_c_grid = comoving_distance(z_grid)；
- geometric_factor = (r_c_z0 / r_c_grid) ** 2；
- alpha = _get_spectral_index() 自动识别幂指数；
- k_correction_factor = ((1 + z0) / (1 + z_grid)) ** alpha；
- total_factor = k_correction_factor * geometric_factor。

随后对每个 z：

如果模型含有红移参数：self._par_z.values = z；
始终按公式修正归一化：self._par_norm.values = original_norm * total_factor[i]。

3.3 谱→计数率（卷积）#

使用 SOXS 完成谱与响应的卷积：
1. spec = soxs.Spectrum.from_pyxspec_model(self._m1)
2. 取能段 0.5–4.0 keV：newspec = spec.new_spec_from_band(0.5, 4.0)
3. 指定响应与曝光：设置 rmf, arf, bkg, exposure, backExposure
4. 用 ARF 卷积得到计数谱：cspec = newspec * soxs.AuxiliaryResponseFile(arf)
5. 总计数率：src_rate = cspec.rate.sum().value
6. 叠加背景贡献（按您的数据口径）：total_rate = src_rate + (bkgnum/duration) * area_ratio
最终计数：total_counts = total_rate * duration

3.4 SNR 计算（Li & Ma公式）#

使用 snr_li_ma(n_src, n_bkg, alpha_area_time)：

n_src = total_counts
n_bkg = bkgnum
alpha_area_time = area_ratio

SNR = sqrt{ 2 [ N_on ln((1+α)/α × N_on/(N_on+N_off)) + N_off ln((1+α) × N_off/(N_on+N_off)) ] }

该公式对高能天文的 on/off 计数统计更稳健。

3.5 自适应搜索 z_max#

在 z_grid 上找到 SNR 首次低于 snr_target 的位置 [z1, z2]；
在 [z1, z2] 区间递归细分，线性插值求解 SNR(z)=目标阈值；
若整个区间 SNR 都高于阈值，逐步扩大上界（max_expand 次）。

4. 推导要点回顾（为何这样做）#

使用共动距离 r_c：

光子数密度是与真实几何距离相关的量，选 r_c 体现 1/r² 衰减；
时间膨胀与能量间隔效应在 XSPEC 光子数谱单位下相互抵消。

K-correction 的来源：

对幂律谱 N ∝ E^-α，能量红移带来额外的 (1+z)^-α 因子；
最终只需考虑该能谱倾斜产生的修正。

统一修正归一化：

无论是否存在 “.Redshift” 参数，归一化都乘以总因子 total_factor，保证谱形在观测能段内的正确缩放。

5. 使用示例#

from autohea.core.utils import RedshiftExtrapolator

ex = RedshiftExtrapolator(
    z0=1.0,
    model="TBabs*zTBabs*powerlaw",
    par=[1e21, 1e21, 1.0, 2.0, 1e-3],  # 示例参数
    arfpath="response.arf",
    rmfpath="response.rmf",
    bkgpath="background.pha",
    srcnum=100,
    bkgnum=1200,
    duration=155,
    area_ratio=1/12,
)

z_max = ex.compute(snr_target=7, show_model_info=True)
print("z_max:", z_max)

6. 常见问题与注意事项#

XSPEC 与 HEASoft：需正确初始化 HEASoft 环境（类内部已调用 HeasoftEnvManager）。
SOXS 版本差异：部分属性在不同版本可能只读，但不影响卷积计算。
模型参数顺序：par 必须与 XSPEC 模型参数顺序完全一致。
红移参数缺失：若模型不含显式红移参数，依旧按 total_factor 修正归一化。
能段选择：示例固定为 0.5–4.0 keV，如需更改可拓展接口。
背景口径：total_rate 中背景项与 area_ratio 的口径请与数据处理保持一致。

7. 验证方法#

使用 verify_redshift_extrapolation(z_test) 可输出详细的几何因子、K 因子与总因子，帮助核对推导：

ex.verify_redshift_extrapolation(z_test=1.5)

打印 r_c、K-correction、总因子与文本解释，便于交叉验证。

8. 后续改进建议#

支持更多谱型的 K-correction 自动识别（如 cutoffpl, grbm 的形状依赖）。
将能段、背景计数处理与响应组合抽象为策略参数，便于复用。
提供误差传播与不确定度评估（如蒙特卡洛抽样）。

如需我们将此文档嵌入到API文档或在Notebook中生成可视化示例，请告知。