Sobol灵敏度指标 — 基于方差分解的全局灵敏度分析

Sobol灵敏度指标的理论基础

Sobol指标的定义

Sobol灵敏度指标是1993年由俄罗斯数学家I.M. Sobol提出的，基于方差分解（ANOVA分解）的全局灵敏度分析方法。与基于局部偏导数的灵敏度分析不同，它评估的是参数空间全域内模型输出的变异性。

ANOVA分解（方差分解）

对于有 $k$ 个独立输入变量 $X_1, X_2, \ldots, X_k$ 的模型输出 $Y = f(X_1, \ldots, X_k)$，ANOVA分解（高维模型表示：HDMR）可写为：

Sobol灵敏度指标的数值计算方法

Saltelli估计量

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}$ 的全部样本估计。

准蒙特卡洛采样

代理模型的结合

• 多项式混沌展开（PCE）：从展开系数能直接分析求出Sobol指标。最高效
• Kriging（高斯过程回归）：预测的不确定性也能定量化，优势明显
• 神经网络：高维、强非线性适配好，但需更多学习数据

Bootstrap置信区间

Sobol指标的估计值有统计不确定性。标准做法是Bootstrap法评估置信区间：

1. 从原 $N$ 样本中有放回地重采样得 $N$ 个新样本
2. 在重采样数据上重新计算Sobol指标
3. 重复 $B$ 次（通常1000次以上）
4. 从得到的 $B$ 个估计值构造百分位置信区间

如果置信区间很宽，需增加样本数 $N$。SALib的 sobol.analyze() 函数中 calc_second_order 选项会自动计算Bootstrap置信区间。

Sobol灵敏度指标的实务应用

Sobol分析的工作流程

1. 问题定义：明确关注的输出指标QoI（Quantity of Interest）和输入参数 $X_1, \ldots, X_k$。确定每个输入的概率分布（均匀分布、正态分布等）
2. 筛选（可选）：参数很多时（20个以上），先用Morris法排除不重要的参数
3. 采样矩阵生成：用Saltelli法生成 $\mathbf{A}, \mathbf{B}, \mathbf{A}_B^{(i)}$。SALib的 saltelli.sample() 最简便
4. 模型评估：在全部采样点运行CAE。推荐并行计算
5. Sobol指标计算：收集评估结果，计算估计量。用SALib的 sobol.analyze()
6. 结果解释和决策：按 $S_i, S_{Ti}$ 排序参数重要性，指导设计优化或不确定性控制

样本规模的确定

结果解释和参数固定

• $S_{Ti} \approx 0$（比如 $< 0.01$） → 该参数可用公称值固定。维数削减
• $S_i$ 很大 → 通过控制该参数的变异（严格公差等）能有效降低输出变异
• $S_{Ti} - S_i$ 很大 → 与其他参数的交互作用强。单独控制效果可能有限

• 板厚: $S_1 = 0.40,\ S_{T1} = 0.48$ — 最支配性
• 屈服应力: $S_2 = 0.25,\ S_{T2} = 0.35$ — 次重要，交互作用也存在
• 摩擦系数: $S_3 = 0.12,\ S_{T3} = 0.18$ — 中等贡献
• 冲压速度: $S_4 = 0.02,\ S_{T4} = 0.03$ — 基本无影响 → 可固定
• 压边力: $S_5 = 0.08,\ S_{T5} = 0.10$ — 小但不能忽视

结论：板厚和屈服应力管理最优先，冲压速度无需考虑，其他要保留。设计指导清晰。

Python SALib实现例

使用SALib（Sensitivity Analysis Library）的最小代码：

import numpy as np
from SALib.sample import saltelli
from SALib.analyze import sobol

# 问题定义
problem = {
    'num_vars': 3,
    'names': ['thickness', 'yield_stress', 'friction'],
    'bounds': [[0.8, 1.2],       # mm
               [200, 350],       # MPa
               [0.05, 0.20]]     # 摩擦系数
}

# Saltelli采样（使用Sobol序列）
param_values = saltelli.sample(problem, 2048, calc_second_order=True)
# → shape: (2048*(3+2), 3) = (10240, 3)

# 模型评估（这里替换为实际CAE求解器）
Y = np.zeros(param_values.shape[0])
for i, X in enumerate(param_values):
    Y[i] = your_cae_model(X[0], X[1], X[2])

# Sobol指标计算
Si = sobol.analyze(problem, Y, calc_second_order=True,
                   conf_level=0.95, print_to_console=True)

# 结果
print("S1:", Si['S1'])       # 第1阶指标
print("ST:", Si['ST'])       # 全阶指标
print("S2:", Si['S2'])       # 二阶交互作用
print("S1_conf:", Si['S1_conf'])  # 95%置信区间

Sobol灵敏度指标的软件对比

主要工具对比

选择指南

• 第一次做Sobol分析 → SALib。pip install就完成，10行代码出结果
• 要高效的PCE方案 → UQLab（MATLAB）或OpenTURNS（Python）
• 大规模HPC任务配合 → Dakota。SLURM/PBS集成充分
• GUI简便操作，Ansys用户 → optiSLang。Workbench直接参数化
• 处理相关输入 → OpenTURNS最灵活

Sobol灵敏度指标的前沿研究

组Sobol指标（Group Sobol Indices）

数学上定义为：对参数组 $u$，$\mathbf{X}_u = \{X_i : i \in u\}$ 的Sobol指标为 $S_u = V[E(Y|\mathbf{X}_u)] / V[Y]$。计算是Saltelli法的扩展。SALib的 groups 参数支持这个。

时间依赖Sobol指标

瞬态分析中，每个时刻 $t$ 都有Sobol指标 $S_i(t), S_{Ti}(t)$。支配参数往往随时间变化，比如碰撞模拟：

• 早期（$t < 5$ ms）：冲击速度 $S_1(t) \approx 0.8$ 支配
• 中期（$5 < t < 20$ ms）：材料参数 $S_2(t) \approx 0.5$ 增加
• 晚期（$t > 20$ ms）：摩擦、边界条件影响增大

处理方法：逐时刻独立做Sobol分析，或用主成分分析（PCA）先降维时序数据再分析。

Shapley效应的对比

Shapley效应 $Sh_i$ 定义为，变量 $X_i$ 在所有可能联盟（子集）中的平均边际方差贡献。保证 $\sum_{i=1}^k Sh_i = V[Y]$ 恒成立，实现"公平分配"。但计算需要 $O(2^k)$ 个子集评估，参数多时爆炸。

Sobol灵敏度指标的故障排除

指标出现负值

1. 样本数 $N$ 不足：估计量是蒙特卡洛估计，$N$ 太小会有统计波动导致负值。增加 $N$（至少加倍）重新计算
2. 该参数真实 $S_i$ 接近零：对输出的贡献很小时，估计噪声会掩盖信号出现负值。看置信区间，如果跨越零，说明"无影响"，不用紧张

$\sum S_i$ 超过1

• 增加 $N$（至少2048以上）
• 确认用的是Sobol序列（准蒙特卡洛）。伪随机收敛慢
• 检查模型输出有没有外点或不连续。仿真中途破裂过

无法收敛

• 输出方差非常大（厚尾分布）→ 试试输出对数变换或标准化
• 输入实际上有相关，但被当成独立处理 → 用Kucherenko法或Shapley效应
• 模型有很强的非线性 → 用代理模型（PCE、Kriging）先平滑再做Sobol
• 模型有随机性噪声（概率模拟）→ 同一输入多次运行求平均

收敛判定的目标：$N$ 加倍时估计值变化 < 10%。