Sobol 灵敏度指数模拟器

Q: 估计 Sobol 指数需要多少次模型评估？

采用 Saltelli 采样时，对于 d 维输入和 N 个基础样本，总计需要 N×(d+2) 次评估即可同时估计 S_i 与 S_Ti（扩展型为 N×(2d+2)）。例如 d=5、N=1024 时约需 7,000 次模型运行。对耗时较长的 CFD/FEM 求解，常见做法是先训练代理模型（克里金、PCE、神经网络），再在代理模型上跑 Sobol。本工具采用的解析公式属于罕见的特例，生产实践中几乎总是使用 Monte Carlo 或拟蒙特卡洛采样。

参数设置

a₁ (X₁ 线性系数)

X₁ 主效应的强度。设为 0 时 X₁ 的主效应消失

a₂ (X₂ 线性系数)

X₂ 主效应的强度

a₃ (X₁X₂ 交互系数)

交互项强度。a₃ = 0 时模型可加，S_T = S 成立

X₁ 的方差 σ₁²

X₁ 的不确定性（假设均值为 0）

X₂ 的方差 σ₂²

X₂ 的不确定性（均值为 0，且与 X₁ 独立）

计算结果

—

总方差 Var[Y]

—

一阶 Sobol S₁

—

一阶 Sobol S₂

—

总 Sobol S_T₁

—

总 Sobol S_T₂

—

交互贡献 S₁₂ (%)

—

ANOVA 分解可视化 — 从输入分布到输出方差

X₁ 与 X₂ 的输入分布经过模型 f(X) 产生输出，其方差按 V₁、V₂、V₁₂ 三部分分解。色块高度表示各 Sobol 贡献的大小。

贡献堆叠条形图 — V₁ / V₂ / V₁₂

扫描曲线 — 总指数 S_T₁ 随交互系数 a₃ 的变化

理论与主要公式

$$Var[Y] = \sum_i V_i + \sum_{i\lt j} V_{ij} + \cdots,\quad S_i = \frac{V_i}{Var[Y]},\quad S_{Ti} = S_i + \sum_{j\neq i} S_{ij} + \cdots$$

ANOVA 分解将总方差拆为主效应 V_i 与交互项 V_{ij}。一阶 S_i 与总 S_Ti 的差值正是输入 i 与其他输入交互作用的重要度量。

$$f(X_1,X_2) = a_1 X_1 + a_2 X_2 + a_3 X_1 X_2,\quad V_1 = a_1^2\sigma_1^2,\; V_2 = a_2^2\sigma_2^2,\; V_{12} = a_3^2\sigma_1^2\sigma_2^2$$

对独立的 X₁、X₂（均值为 0），各分量可写为闭式。当 a₃ = 0 时模型变为可加，S_T 与 S 相等。

Sobol 指数与全局灵敏度分析

🙋

灵敏度分析不就是稍微改一下输入看输出怎么变吗？Sobol 指数有什么特别的呢？

🎓

问得好。本科里学的灵敏度大多是 ∂f/∂x 那种「在某一点处的局部灵敏度」。Sobol 属于另一类，叫「全局灵敏度分析（GSA）」。它把每个输入 X_i 当作随机变量，问：在整个输入分布上，输出 Y 的方差有多少是由 X_i 引起的？答案是 0～1 之间的贡献率，即使 f 高度非线性、非可加也仍然有意义。这就是 CFD 和气候代码把它当作标配工具的原因。

🙋

默认设置下卡片上写着 S₁=0.643 和 S_T₁=0.714。这两个值为什么不一样？

🎓

这正是整个面板里最有信息量的差异。S_i 表示「把 X_i 固定为某个值时，输出方差能减少多少」。S_Ti 表示「除 X_i 外都固定时，剩下多少方差」。差值 S_T₁ − S₁ = 0.071 恰好是「X₁ 通过与 X₂ 的交互而贡献的方差份额」。把 a₃ 滑到 0 试试，你会看到 S₁ 与 S_T₁ 重合，这就是可加模型的标志。

🙋

那么把交互系数 a₃ 调大，X₁ 通过交互项「借来」更多影响力，S_T₁ 就跟着变大对吧？

🎓

完全正确。下方的扫描曲线就显示 S_T₁ 会随 |a₃| 增大而趋近 1。物理上的图景是「X₁ 的效应取决于 X₂ 当前的取值」。经典例子是机翼升力系数对迎角和马赫数的依赖——∂C_L/∂α 本身就随 M 变化。只看 S₁ 会说「迎角中等重要」，而看 S_T₁ 才能揭示它的真实重要度。

🙋

真实的 CFD 或 FEM 模型可没有这么干净的解析公式。实际工程里大家怎么算 Sobol？

🎓

靠采样。标准估计器是 Saltelli 法：d 维输入、基础样本数 N 的情况下，需要约 N·(d+2) 次模型评估即可同时得到 S_i 和 S_Ti。比如 d=10、N=1024，大约要 12,000 次评估。一次 CFD 算一个小时的话，那就是 500 天起步。所以生产流程一般先训练代理模型——克里金、多项式混沌展开（PCE）或神经网络——再在便宜的代理上跑 Sobol。先用 Morris 做筛选、再用 Sobol 精算的两阶段法是最具性价比的标准做法。

🙋

如果某个输入的 S_T 非常小，直接把它从模型里删掉就行吗？

🎓

大多数情况下可以——这就是 Sobol 最经典的「降维」用途。S_Ti ≈ 0 的输入可以固定在名义值，从后续研究中剔除。但要注意两点：第一，这个结论只在「当前输入分布范围」内成立，把范围扩大后原本很小的 S_T 可能会显著上升；第二，Sobol 是基于方差的，对那种「只改变均值、不改变离散度」的输入会视而不见。涉及安全系数、设计许用值等关注均值的量时，务必再用矩独立指标或 Borgonovo δ 复核一次。盲目按「S_T 小所以可以扔」操作的团队，吃过亏的不在少数。

常见问题

一阶指数 S_i 表示输入 X_i 单独作用（主效应）对输出方差的贡献率，S_i = V_i/Var[Y]。总指数 S_Ti 包含所有涉及 X_i 的项（主效应＋全部交互），即 S_Ti = S_i + Σ_{j≠i} S_{ij} + …。两者之差 S_Ti − S_i 反映 X_i 与其他输入的交互强度，差越大说明仅固定 X_i 一项无法把输出方差降低到主效应所暗示的程度。

Morris 基本效应法以及偏相关系数 PCC、标准化回归系数 SRC 计算便宜，但只看输入分布上的局部灵敏度。当输出非线性、非可加时，这些局部指标可能给出错误排序。Sobol 属于全局灵敏度分析（GSA），在整个输入概率分布上积分，不假设可加性或线性。代价是收敛需要约 N·(d+2) 次评估（Saltelli 估计），成本较高。务实做法是：先用 Morris 做筛选，再对剩下的少量重要变量执行 Sobol。

采用 Saltelli 采样时，对于 d 维输入和 N 个基础样本，总计需要 N·(d+2) 次评估即可同时估计 S_i 与 S_Ti（扩展型为 N·(2d+2)）。例如 d=5、N=1024 时约需 7,000 次模型运行。对耗时较长的 CFD/FEM 求解，常见做法是先训练代理模型（克里金、PCE、神经网络），再在代理模型上跑 Sobol。本工具采用的解析公式属于罕见的特例，生产实践中几乎总是使用 Monte Carlo 或拟蒙特卡洛采样。

若 S_Ti 接近 0，则该输入对输出方差几乎没有贡献，在 UQ 范围内可以固定为名义值，这就是 Sobol 经典的筛选与降维用途。但有两点需要注意：第一，该结论仅在当前输入分布范围内成立，扩大分布范围后小的 S_T 可能急剧上升；第二，S_T 是基于方差的指标，会忽略仅影响输出均值（偏置）而不改变离散度的输入。对安全系数、设计许用值等关注均值的量，应结合矩独立指标或 Borgonovo δ 等额外指标确认。

现实应用

CFD 设计优化：将翼型阻力系数 C_D 作为输出，把迎角、马赫数、雷诺数、湍流模型常数等 10～30 个不确定输入交给 Sobol，可以揭示真正主导设计的参数。例如「马赫数的 S_T 高达 0.7」这样定量的排序，可直接用于缩小设计空间、给代理模型采样排定优先级。常见组合是 OpenFOAM + Dakota、SU2 + UQpy。

FEM 结构许用设计：材料常数（杨氏模量、屈服应力、密度）与几何尺寸（壁厚、圆角半径）都有公差。Sobol 把它们对峰值应力或固有频率方差的贡献分解出来，工程师可对少数主导输入加严公差、对其余输入放宽至规格上限以降低成本。航空与汽车的认证文件越来越多地要求附上这种 UQ 摘要表。

气候与环境模型：大气环流模式 GCM 含有数百个参数化常数，Sobol 指数是溯源云量或平衡气候敏感度（ECS）不确定性的主力工具。一次 GCM 运行动辄几小时到数天，故 Sobol 几乎总是在 PCE 代理上完成。IPCC 工作组报告中的多项不确定性结论正是基于这种方差分解。

可解释机器学习：近期研究使用 Sobol 作为 SHAP 的补充，对神经网络或梯度提升模型的特征重要度进行排序。SHAP 擅长局部贡献，而 Sobol 给出的是在整个输入分布上的贡献，更适合记录模型在受监管领域（生物医药、医疗、金融）中的整体行为。

常见误解与注意点

最大的陷阱是「在相关输入下直接套用 Sobol」。经典估计器假设输入相互独立，但实际工程里「降水量和气温」「屈服应力与杨氏模量」往往是相关的。强行使用独立假设会让相关输入互相蚕食贡献，得到毫无意义的排序。对策是改用相关版扩展（Kucherenko 方法、Mara & Tarantola 的 Shapley 效应），或先用 PCA 将输入变换为无关的主成分再做灵敏度分析。看到一份没有声明独立性如何验证的 Sobol 表格，请保持怀疑。

其次是「不检查采样收敛性」。Saltelli 估计仅在 N→∞ 时收敛到真值，有限 N 下存在偏差，尤其 S_T 偏差最为顽固。N=128 时算出 S_T₁=0.3，N=4096 时变成 S_T₁=0.6 的情况屡见不鲜。可信的研究会把 N 翻倍直到 S_i 与 S_Ti 稳定，并用 bootstrap 给出置信区间。若一篇论文只把 Sobol 数字摆上幻灯片而没有收敛图和置信区间，基本可以不必当真。

最后，「把 Sobol 当作唯一指标」也是常见错误。Sobol 仅基于方差，描述输出分布的「离散度」而非完整形状。对厚尾分布（绝大多数时刻安全，极少数情形会发生灾难性后果），Sobol 排名靠后的输入反而可能主导尾部风险。现代 UQ 通常把 Sobol 与矩独立指标（Borgonovo δ）、PAWN 密度法、VaR/CVaR 等尾部分位灵敏度结合使用。Sobol 强大但并非全能，要把它视为灵敏度仪表盘上的一项指标，而非唯一指标。

Sobol 指数与全局灵敏度分析

常见问题

现实应用

常见误解与注意点

使用指南

具体计算示例

实务注意事项