CF-LIBS详细步骤

介绍实现CF-LIBS的详细步骤

本文定位：记录具体实现中的测温、密度估计、配分函数、Saha 合并、谱线筛选与后处理。

关联阅读： 基本原理与流程 · 物理假设与诊断 · 指标解释与结果表达。

CF-LIBS 执行流程详细拆解

CF-LIBS（Calibration-Free LIBS，GB/T 38257—2019 第 9.2.3 条）的核心逻辑位于 libs_analysis/src/libs_analyzer/models/cf_libs.py 的 CFLIBS 类中。整个多元素全成分输出由 predict_composition() 编排，依次执行以下 5 个核心步骤 + 1 个后处理步骤。

Step 1: Boltzmann 图测温

方法: _estimate_temperature(X)

物理原理: LTE 假设下，同一元素多条谱线强度遵循 Boltzmann 分布。对目标元素（默认 Fe I）的各谱线取对数：

以 $E_j$ 为自变量做线性回归，斜率 $= -1/(k_B T) \Rightarrow T$。

代码流程:

1. 加载测温谱线: 取 temp_element（默认 Fe）temp_ion（默认 I）的谱线表；若 temp_min_wl_nm > 0（默认 400nm），剔除紫外区谱线（避免 Fe 密集谱区宽窗背景被邻线抬高导致 T 虚高）。
2. 要求 ≥ max(min_lines_for_boltzmann, 4) 条线: 4 条是统计稳定的下限（2 点拟合 R² 恒为 1，门控失效）。不足则全部样品回退 t_ref_k。
3. 调用 _multi_line_net_and_snr() 提取净强度/SNR/对比度:
   • 净强度 net: 窗内峰值（或面积）减去宽窗中位数背景。
   • 噪声: 背景区 IQR / 1.349（对 clip(≥0) 退化分布鲁棒，不用 MAD）。
   • SNR: (peak − bg) / noise。
   • 对比度: peak / shoulder_median，shoulder 为紧邻峰窗两侧的局部中位强度，用于过滤连续背景上的伪峰。
4. 计算 $y = \ln(\text{net} / gA)$，逐样本线性回归。
5. 谱线过滤:
   • 主过滤：SNR ≥ snr_threshold（默认 5.0）且对比度 ≥ min_peak_contrast（默认 2.0）。
   • 放宽：不足时退化为仅 SNR ≥ 5.0；仍不足则回退 T_ref。
6. 强线自吸收规避: 取净强度 ≤ temp_line_quantile（默认 0.3，即最弱 30%）分位数的弱线子集。物理依据：线光学深 τ ∝ 线强，强线被自吸收压低 → Boltzmann 斜率变浅 → T 虚高 5000–10000K。弱线不足 4 条时回退全部有效线。
7. 能量跨度预检: $\Delta E = \max(E) - \min(E) < $ min_boltzmann_delta_e_ev（默认 1.2 eV）时回退（斜率对误差过度敏感）。
8. 初拟合: _boltzmann_fit(E, y) → (slope, R², intercept)。
9. 向下离群点移除（自吸收）: 迭代移除残差最低的点（残差 $< -$ self_absorb_sigma × RMSE），最多移除 30% 且移除后 ≥ 4 条。仅移除向下偏离（自吸收只压低强度），不删向上离群点（避免旧贪心算法过拟合到 2 点拟合的退化问题）。
10. 移除端点后复检 ΔE（可能缩小）。
11. 斜率方向检查: slope 需为负（$-1/kT < 0$）。
12. R² 显著性门控: 用 F 检验 95% 临界值动态调整 R² 下限（n=4 需 R²≥0.90，n=5 需 0.77，n=6 需 0.66），取 max(min_boltzmann_r2, r2_sig)。scipy 缺失时退化为固定门槛。
13. 温度合理性: $T \in [$ t_min_k, t_boltzmann_max_k $]$（默认 [3000, 20000] K）。此判定在删线结束后才执行，不参与删线决策，避免"删线直到温度过线"的过拟合。
14. 不满足任一条件的样品回退 t_ref_k（默认 10000 K）。
15. 最终裁剪到 [t_min_k, t_max_k]。

Step 2: 电子密度估计（Stark 展宽法）

方法: _estimate_electron_density(X, T) → 调用 electron_density.estimate_ne()

物理原理: 等离子体中带电粒子间的电场导致谱线 Stark 展宽，Lorentzian 分量 FWHM 与电子密度成正比：

其中 $w$ 为 Stark 展宽参数（nm）。

代码流程:

1. 优先级判定:
   • electron_density 为 float → 用户指定固定值（最高优先，不受 enable_stark_density 开关影响）。
   • enable_stark_density=False（消融实验）→ 固定 DEFAULT_NE_CM3 = 1×10^{17} cm⁻³。
   • 否则 → Stark 展宽法自动估计。
2. Stark 展宽法自动估计 (estimate_ne()):
   • 加载 Stark 表: 从 YAML 读取候选谱线的波长、Stark 参数 $w$、共振标志、多重线信息。
   • Voigt 拟合: 对每条候选线、每个样品，在 ±fit_window_nm（默认 1.0nm）窗内拟合 Voigt + 线性背景（_voigt_bg()），提取 Lorentzian HWHM $\gamma$。
   • 仪器展宽扣除: 用 Olivero-Longbothom 1977 近似反卷积，从观测 Voigt FWHM 中扣除仪器高斯展宽（instrument_width_nm，默认 0.08nm），得到纯 Stark FWHM。
   • 逐线计算 $N_e$: $N_e = 10^{17} \times \text{FWHM}_{\text{Stark}} / w$。
   • 自吸收检测: 对共振线（如 Mg II 279/280、Ca II 393/396）检查多重线强度比偏离理论值的程度，超过 self_absorption_tolerance（默认 0.4）的样品标记为自吸收，该线不参与 $N_e$ 计算。
   • 中位数聚合: 非共振线中位数优先 → 通过自吸收检查的共振线中位数 → 回退 1×10^{17}。
   • 物理合理性裁剪: $N_e \in [10^{15}, 10^{19}]$ cm⁻³。
3. 失败回退: 异常时回退 1×10^{17} cm⁻³。
4. 返回: (Ne (n,), diagnostics: {ne_source, per_line, n_lines_used, ...})。

Step 3: 配分函数计算

方法: _Z(element, ion, T) → 调用 partition_function.partition_function()

物理原理: 配分函数 $Z_{s,\text{ion}}(T)$ 是元素在不同能级上布居的统计权重之和，是 CF-LIBS 浓度公式中消去 $F$ 后唯一与温度相关的修正因子。不同元素、不同电离态的 $Z(T)$ 差异显著（如 Fe I 在 10000K ≈ 28，Fe II ≈ 42）。

代码流程:

1. 优先级判定:
   • 用户提供了 partition_functions（旧格式 {Fe: 8.0} 固定值）且未指定 partition_functions_path → 直接返回固定值，忽略 ion 与 T（向后兼容）。
   • 否则 → 调用 NIST 表插值。
2. NIST 表插值 (partition_function()):
   • 加载表: load_partition_table() 从 configs/partition_functions.yaml 读取 {T_points: [...], elements: {elem: {ion: {T: Z}}}}，校验 Z ∈ [0.5, 1000] 且随温度单调递增。
   • 线性插值: 对输入 T（可为数组）在 NIST 离散温度点上做线性插值。T 超出表范围时用边界值外推。
   • 缺失回退: 元素/电离态不在表中时返回 fallback=1.0。
3. 返回: 与 T 同 shape 的 ndarray。

Step 4: Saha 电离平衡合并

方法: 在 _compute_relative_concentrations() 内部执行 → 调用 saha.saha_ratio()

物理原理: Saha 方程描述相邻电离态的浓度比：

CF-LIBS 浓度公式对每个电离态分别给出 $q_{\text{ion}} = F \cdot C_{s,\text{ion}} / Z_{s,\text{ion}}(T)$，需合并得到元素总浓度 $C_s = C_{s,\text{I}} + C_{s,\text{II}}$。

代码流程（在 _compute_relative_concentrations() 中对每个元素执行）:

1. 枚举可用电离态: _available_ions(elem) 返回原子参数表中存在的 I/II 态（跳过 III 及以上，LIBS 温度下占比 <1%）。
2. 各电离态分别计算 $q_{\text{ion}}$: 对每个 ion，取谱线表 → 提取 net/SNR/contrast → 多重过滤 → 计算 $q = \text{net} / (gA \cdot \exp(-E/k_BT))$。
3. 选择参考电离态: _pick_reference_ion(ions) 优先选 I（中性原子谱线更全），单态则直接用。
4. 合并策略（enable_saha_merge=True 时）:
   • 两态均有谱线: $C_{\text{total}} = C_{\text{ref}} + C_{\text{other}}$（直接相加，$F$ 相同在归一化时消去）。
   • 仅参考态有谱线: 用 Saha 方程推算另一态。
     • 参考态为 I → $C_{\text{total}} = C_{\text{I}} \times (1 + \text{saha\_ratio})$
     • 参考态为 II → $C_{\text{total}} = C_{\text{II}} \times (1 + 1/\text{saha\_ratio})$
   • 放大钳位: saha_max_amplification（默认 $10^4$）。T 过高/Ne 过低时放大倍数可达 $10^3$–$10^4$，超限时钳位并记录诊断（saha_clipped）。
5. 单电离态或 Saha 关闭: 直接 $C_{\text{raw}} = q \times Z(T)$，跳过合并。

Step 5: 多重谱线过滤 + 相对浓度计算 + 归一化

方法: _compute_relative_concentrations(X, T, Ne)

这是整个流程中最复杂的步骤，包含谱线质量筛选、浓度计算、降级检测和归一化。

物理原理:

归一化消去总比例因子 $F$：$x_s = C_{s,\text{raw}} / \sum_{s'} C_{s',\text{raw}}$。

代码流程:

5a. 逐元素、逐电离态的谱线质量过滤

对每个元素的每个电离态，执行以下 五层过滤（每层可由 enable_* 开关独立关闭，用于消融实验）：

5b. 浓度 $q$ 计算（三条路径）

对每个样品 $i$，统计通过所有过滤的有效线数 $n_{\text{valid}}$：

5c. Saha 合并（见 Step 4）

5d. 降级元素浓度上限

对走了降级路径的元素（且未被任一电离态主路径覆盖的样品），施加地质合理浓度上限（trace_concentration_caps，如 Mn ≤ 0.5%、S ≤ 0.8% 等），防止噪声/干扰放大致过估。仅约束降级元素，主量元素不受影响。 钳位后重新归一化保持总和为 1。

5e. 归一化

$x_s = C_{s,\text{raw}} / \sum_{s'} C_{s',\text{raw}}$，消去未知总比例因子 $F$，得到各元素原子摩尔分数。

Step 6（后处理）: LTE 校验 + 原子浓度转氧化物质量分数

方法: mcwhirter_check() + atomic_to_composition()

6a. McWhirter LTE 判据校验

调用: mcwhirter_check(Ne, T, delta_E)

物理原理: LTE 假设成立需满足 McWhirter 判据：

其中 $\Delta E$ 取参与浓度计算的谱线中最大上能级 $E_{\text{upper}}$（保守上界，因跃迁能量 ≤ $E_{\text{upper}}$）。

代码流程:

• 计算各样品的阈值并与实测 $N_e$ 比较。
• 返回 {satisfied: (n,) bool, threshold: (n,), min_ne: float, delta_E_ev: float}。
• 不满足的样品不代表结果一定错误，但置信度降低（国标要求报告此判据）。

6b. 原子摩尔分数 → 氧化物质量分数

调用: atomic_to_composition(C_atomic, feo_fe_ratio)

物理原理: CF-LIBS 输出原子摩尔分数，但铁矿石分析报告要求氧化物/元素质量分数。转换需：

1. 摩尔分数 × 原子量 → 元素质量。
2. 按化学计量比配氧（如 Si → SiO₂ 需 2 个 O 原子）。
3. 归一化到全氧化样品总质量（金属 + 氧 = 100%）。
4. 乘以氧化物转换因子（如 Si → SiO₂ 的 $M_{\text{SiO}_2}/M_{\text{Si}} = 2.1392$）。
5. P 和 S 以元素形式报告（铁矿石国标惯例）。
6. FeO 估算：$\text{FeO} = \text{Fe}^{2+} \times M_{\text{FeO}} / M_{\text{Fe}}$，其中 $\text{Fe}^{2+} = \text{Fe}_{\text{total}} \times $ feo_fe_ratio（磁铁矿 ≈ 0.333，赤铁矿 ≈ 0）。

代码流程:

1. atomic_to_element_mass(): 原子摩尔分数 → 含氧元素质量分数（氧按化学计量比分配到各金属元素的分母中）。
2. element_mass_to_composition(): 元素质量分数 → 10 组分 wt% 报告（TFe、FeO、SiO₂、Al₂O₃、CaO、MgO、TiO₂、MnO、P、S）。

总流程总结

predict_composition(X)
  │
  ├── Step 1: _estimate_temperature(X)
  │     Boltzmann 图法逐样品测温 → T (n,)
  │     核心机制：弱线选择规避自吸收 + 向下离群点移除 + F 检验 R² 门控
  │
  ├── Step 2: _estimate_electron_density(X, T)
  │     Stark 展宽法估计 → Ne (n,)
  │     核心机制：Voigt 拟合 + 仪器展宽扣除 + 共振线自吸收检测 + 中位数聚合
  │
  ├── Step 3: _Z(element, ion, T)  [在 Step 5 内逐元素调用]
  │     NIST 表插值 → Z(T) (n,)
  │     核心机制：离散温度点线性插值，缺失回退 1.0
  │
  ├── Step 4: Saha 合并  [在 Step 5 内逐元素执行]
  │     saha_ratio() → C_II/C_I → C_total
  │     核心机制：两态直接相加 or 单态 Saha 推算 + 放大钳位
  │
  ├── Step 5: _compute_relative_concentrations(X, T, Ne)
  │     五层过滤 + 三条浓度路径 + Saha 合并 + 降级上限 + 归一化
  │     → {element: 原子摩尔分数 (n,)}
  │
  └── Step 6: 后处理
        ├── mcwhirter_check(Ne, T, ΔE)  → LTE 判据
        └── atomic_to_composition(C_atomic) → 10 组分 wt% 报告

calibrated 模式（可选）

若 calibrated=True，在 fit(X, y) 阶段用已知浓度 $y$ 对目标元素（默认 Fe）的相对浓度 $c_{\text{rel}}$ 做线性校正 $y = a \cdot c_{\text{rel}} + b$。predict_composition() 中以该校正比等比缩放全部元素后重新归一化——即国标 9.2.3 允许的"单点/多点定标 CF-LIBS"模式。