Sobol灵敏度指标——基于方差分解的全局灵敏度分析
理论与物理
Sobol指标是什么
老师，Sobol灵敏度指标能让我们了解什么？灵敏度分析有很多种类吧。
简单来说，它是一种量化输入参数对输出变异的贡献程度的方法。一阶Sobol指标 $S_i$ 表示主效应，全阶Sobol指标 $S_{Ti}$ 表示包含交互作用在内的贡献率。
“贡献率”具体是什么意思呢？
例如，如果 $S_i = 0.8$，就意味着仅该参数就能解释输出方差的80%。即使有10个设计参数，真正重要的往往只有2~3个。Sobol指标就是用数值来展示这一点。
诶，有10个设计变量但重要的只有2~3个吗！如果知道了这一点，剩下的7~8个就可以固定了吧。
正是如此。在CAE实际工作中，比如汽车碰撞模拟中有板厚、材料特性、摩擦系数、焊接条件……等大量输入参数，但如果通过Sobol指标知道“仅板厚和材料的屈服应力就能解释输出变异的90%”，那么其他参数就可以用标称值固定。计算成本会大幅减少。
Sobol灵敏度指标是1993年由俄罗斯数学家 I.M. Sobol 提出的，基于方差分解（ANOVA分解）的全局灵敏度分析方法。与基于局部偏微分的灵敏度分析不同，它评估的是参数空间全域内模型响应的变异。
ANOVA分解（方差分解）
对于具有 $k$ 个独立输入变量 $X_1, X_2, \ldots, X_k$ 的模型输出 $Y = f(X_1, \ldots, X_k)$，ANOVA分解（高维模型表示: HDMR）可以写成：
其中 $f_0 = E[Y]$ 是总均值，$f_i(X_i)$ 是变量 $X_i$ 的主效应，$f_{ij}$ 是二阶交互作用项。为了使此分解唯一，需要各项相互正交（积分为零）的条件，前提是输入变量相互独立。
方差也可以类似地分解：
每个 $V_i, V_{ij}, \ldots$ 是对应项的方差，表示对总方差 $V[Y]$ 的部分方差贡献。
一阶Sobol指标 $S_i$
具体是什么数学公式呢？
一阶Sobol指标是，固定 $X_i$ 时的条件期望的方差，除以输出总方差。
其中，
- $E_{X_{\sim i}}(Y \mid X_i)$：固定 $X_i$，对所有其他变量 $X_{\sim i}$ 取平均后的条件期望
- $V_{X_i}[\cdot]$：该条件期望关于 $X_i$ 的方差
- $V[Y]$：输出 $Y$ 的无条件方差（总方差）
$S_i$ 的取值范围是 0~1，$S_i = 0$ 表示 $X_i$ 对输出完全没有影响。$\sum_{i=1}^{k} S_i = 1$ 意味着交互作用为零（可加模型）。
“对所有其他变量 $X_{\sim i}$ 取平均”是关键所在呢。意思是只变动 $X_i$，而将其他所有变量平均化，看输出有多大波动，对吧。
理解完全正确。这是从总方差定理 $V[Y] = V_{X_i}[E_{X_{\sim i}}(Y|X_i)] + E_{X_i}[V_{X_{\sim i}}(Y|X_i)]$ 推导出来的公式。第一项是“可由 $X_i$ 的主效应解释的方差”，第二项是“即使固定 $X_i$ 后仍剩余的方差”。
全阶Sobol指标 $S_{Ti}$
$S_{Ti}$ 表示 $X_i$ 参与的所有交互作用项在内的综合贡献率。一阶指标与全阶指标的差值 $S_{Ti} - S_i$ 越大，说明 $X_i$ 与其他变量的交互作用越强。
能用一个实际例子说明 $S_i$ 和 $S_{Ti}$ 的区别吗？
例如，某个换热器设计有3个输入（流速 $X_1$、管径 $X_2$、材料导热率 $X_3$）。结果是：
- $S_1 = 0.45,\ S_{T1} = 0.52$ — 流速主效应大，交互作用较小
- $S_2 = 0.10,\ S_{T2} = 0.35$ — 管径主效应小，但与其他变量的交互作用大
- $S_3 = 0.30,\ S_{T3} = 0.33$ — 导热率几乎只有主效应
管径 $X_2$ 如果只看 $S_2$ 容易判断为“影响小”，但看 $S_{T2}$ 则发现不可忽视。流速和管径的组合（类似雷诺数效应）以交互作用的形式体现出来了。
原来如此！只看 $S_i$ 就断定“这个变量不需要”的话，可能会有危险的情况呢。
所以，判断参数是否可固定时，必须使用 $S_{Ti}$。如果 $S_{Ti} \approx 0$ 就可以安全地固定。不能仅凭 $S_i$ 判断。这是实际工作中非常重要的要点。
高阶交互作用指标
二阶交互作用指标定义如下：
其中 $S_{ij}^{\text{closed}} = \frac{V_{X_i, X_j}[E_{X_{\sim ij}}(Y|X_i, X_j)]}{V[Y]}$ 是 $X_i$ 和 $X_j$ 的闭合二阶指标。理论上可以计算任意阶次，但由于计算成本呈指数增长，实际工作中通常 $S_i$ 和 $S_{Ti}$ 两个指标就足够了。
Sobol指标性质总结
- $0 \leq S_i \leq S_{Ti} \leq 1$：一阶指标小于等于全阶指标
- $\sum_{i=1}^k S_i \leq 1$：等号在交互作用为零（可加模型）时成立
- $\sum_{i=1}^k S_{Ti} \geq 1$：等号在交互作用为零时成立
- 以输入独立性为前提：对于相关的输入，需使用广义Sobol指标（Kucherenko 2012）或Shapley效应
与局部灵敏度分析的区别
- 局部灵敏度分析（偏微分 $\partial Y/\partial X_i$）：仅捕捉基准点附近的线性响应。计算成本低，但无法看到非线性效应和交互作用
- 全局灵敏度分析（Sobol）：评估参数空间全域的变异。能捕捉所有非线性和交互作用。但需要蒙特卡洛式采样，计算成本高
- Morris法（筛选法）：适用于“重要或不重要”的二值判断。常作为Sobol指标的前置步骤使用
数值解法与实现
Saltelli估计量
定义式我明白了，但实际怎么计算 $S_i$ 和 $S_{Ti}$ 呢？那个条件期望的方差，解析地求不出来吧？
问得好。实际工作中使用Saltelli (2002, 2010) 提出的蒙特卡洛估计量。首先生成两个独立的样本矩阵 $\mathbf{A}$ 和 $\mathbf{B}$（各 $N \times k$）。然后，构造一个矩阵 $\mathbf{A}_B^{(i)}$——这是将 $\mathbf{A}$ 的第 $i$ 列替换为 $\mathbf{B}$ 的第 $i$ 列得到的。
Saltelli估计量的步骤：
步骤 1：独立生成 $N$ 行 $k$ 列的样本矩阵 $\mathbf{A}$ 和 $\mathbf{B}$（推荐使用准蒙特卡洛序列）。
步骤 2：对每个 $i = 1, \ldots, k$，构造混合矩阵 $\mathbf{A}_B^{(i)}$。第 $j$ 列为：
步骤 3：对 $\mathbf{A}$, $\mathbf{B}$, $\mathbf{A}_B^{(1)}, \ldots, \mathbf{A}_B^{(k)}$ 的所有行评估模型 $f$。总共需要 $N(k+2)$ 次模型评估。
步骤 4：计算估计量。
其中 $V[Y]$ 是从 $\mathbf{A}$ 和 $\mathbf{B}$ 的全部样本中估计的。
也就是说，生成 $\mathbf{A}$ 和 $\mathbf{B}$ 两套样本，然后一列一列地替换，来调查“只改变这个变量会怎么样”，对吧。
对，正是如此。多亏了这个“列替换”技巧，才能用蒙特卡洛高效地估计定义式中的多重积分。顺便说一下，对于 $S_{Ti}$ 的估计式，还有Jansen (1999) 或Sobol-Jansen估计量等变体，据说Jansen估计量的偏差更小。上面写的是Jansen估计量。
准蒙特卡洛采样
在生成样本矩阵那里写着“推荐准蒙特卡洛序列”，用普通的随机数不行吗？
可以用但收敛慢。伪随机数的蒙特卡洛收敛速度是 $O(1/\sqrt{N})$，而使用Sobol序列或Halton序列等准蒙特卡洛法可以改善到 $O((\log N)^k / N)$。实际工作中，通常可以用1/10左右的样本数达到相同精度。
| 采样方法 | 收敛速度 | 特点 | 推荐场景 |
|---|---|---|---|
| 伪随机数MC | $O(N^{-1/2})$ | 实现简单，存在聚类 | 维数少、计算轻量时 |
| Sobol序列（QMC） | $O(N^{-1}(\log N)^k)$ | 低差异度，高维也有效 | 标准推荐。SALib默认 |
| Halton序列 | 同上 | 高维时质量易劣化 | 低维（$k \leq 10$）适用 |
| 拉丁超立方（LHS） | $O(N^{-1/2})$程度 | 边缘分布采样准确 | Sobol序列无法使用时的替代 |
与代理模型并用
如果一次CAE模拟需要几个小时，那么评估 $N(k+2)$ 次模型不是很困难吗？参数10个，$N=1000$ 的话就是12,000次…
这正是实际工作中的瓶颈。解决方案是使用代理模型。用较少的CAE评估点（数十到数百点）构建响应面，然后将大量样本输入代理模型来估计Sobol指标。
- 多项式混沌展开（PCE）：可以从展开系数直接解析计算Sobol指标。最有效率
- Kriging（高斯过程回归）：优势在于也能量化预测的不确定性
- 神经网络：适用于高维、强非线性，但需要大量训练数据
PCE可以直接从展开系数得到Sobol指标吗！那很方便啊。
在PCE展开式 $Y \approx \sum_{\boldsymbol{\alpha}} c_{\boldsymbol{\alpha}} \Psi_{\boldsymbol{\alpha}}(\mathbf{X})$ 中，一阶指标 $S_i = \sum_{\boldsymbol{\alpha} \in \mathcal{A}_i} c_{\boldsymbol{\alpha}}^2 / \sum_{\boldsymbol{\alpha} \neq \mathbf{0}} c_{\boldsymbol{\alpha}}^2$ 表示为系数平方和的比值。这样就无需蒙特卡洛采样，速度极快。UQLab 和 OpenTURNS 支持这种方法。
Bootstrap置信区间
Sobol指标的估计值伴随有统计不确定性。使用Bootstrap法评估置信区间是标准做法：
- 从原始 $N$ 个样本中有放回地重采样 $N$ 个
- 用重采样数据重新计算Sobol指标
- 重复此过程 $B$ 次（通常1000次以上）
- 根据得到的 $B$ 个估计值构建百分位置信区间
如果置信区间较宽，则需要增加样本数 $N$。在SALib中，通过 sobol.analyze() 的 calc_second_order 选项会自动计算Bootstrap置信区间。
实践指南
Sobol分析的工作流程
实际进行Sobol分析时，从头到尾需要做什么呢？
6个步骤。
- 问题定义：明确输出QoI（关注量）和输入参数 $X_1, \ldots, X_k$。确定每个输入的概率分布（均匀分布、正态分布等）
- 筛选（可选）：参数数量多时（20个以上），先用Morris法排除不重要的参数
- 生成样本矩阵：用Saltelli法生成 $\mathbf{A}, \mathbf{B}, \mathbf{A}_B^{(i)}$。使用SALib的
saltelli.sample()最方便 - 模型评估：在所有生成的样本点上运行CAE模型。推荐并行计算
- 计算Sobol指标：收集评估结果并计算估计量。使用SALib的
sobol.analyze() - 结果解释与决策：根据 $S_i, S_{Ti}$ 对参数重要性进行排序，并反映到设计优化或不确定性降低中
样本大小的确定方法
$N$ 应该取多少呢？越大越好我明白，但预算也有限…
我总结一个参考表吧。
| 参数数量 $k$ | $N$（基本样本数） | 总模型评估次数 $N(k+2)$ | 备注 |
|---|---|---|---|
| 3〜5 | 1,024〜2,048 | 5,120〜14,336 | 直接Sobol法足以应对 |
| 5〜10 | 2,048〜4,096 | 14,336〜49,152 | 接近推荐使用代理模型的界限 |
| 10〜20 | 4,096〜8,192 | 49,152〜180,224 | 必须使用代理模型。推荐事先进行Morris筛选 |
| 20+ | — | — | 直接Sobol法不现实。采用PCE或Morris→Sobol的两阶段策略 |
$N$ 最好是2的幂次吗？
使用Sobol序列（准蒙特卡洛）时，$N = 2^p$ 是最优的。可以最大限度地利用其低差异度的性质。如果是伪随机数，任意 $N$ 都可以。
结果的解释与参数固定
得到Sobol指标后，具体如何进行判断呢？
有三个判断标准：
- $S_{Ti} \
相关主题
この記事の評価ご回答ありがとうございます！参考に
なったもっと
詳しく誤りを
報告