2026/9/8 0:35:20
FEATURED · 精选文章 
做敏感性分析最怕的不是算不出来而是算了一堆数最后没法回答“这模型到底该调哪个参数”。我之前在做一个带七八个可调参数的仿真模型时跑一遍完整工况要一两分钟手动调参根本试不过来做局部分析又发现同一个参数在不同基准点上影响完全不一样——这时候才意识到需要的不是某个点附近的导数而是参数在整个取值空间里对输出的影响排序。找了一圈之后锁定了eFAST扩展傅里叶幅度敏感性检验。这个方法的优势在于给每个参数分配一个“扫描频率”把参数空间转换成周期信号来考察输出方差——既能给出每个参数单独的影响一阶指数也能给出包含交互作用的总影响总指数而且不像Sobol那样动辄需要几万次模型运行。不过MATLAB里没有现成的eFAST函数命令行敲efast是什么干货都搜不到的。我花了一周左右把原理吃透、代码写完、再用标准测试函数验证了一遍整个过程踩了不少坑。这篇就按我实际操作的顺序把方法选型、核心原理、MATLAB实现代码、验证结果、以及调参避坑经验完整写出来给同样需要用全局敏感性分析的朋友一条能直接照做的路径。1. 为什么是eFAST先厘清敏感性分析的几种路线1.1 局部敏感性分析与全局敏感性分析的取舍很多人在模型调参初期都会试着做一遍“控制变量法”固定其他参数把某个参数从最小值调到最大值看输出怎么变。这就是最典型的局部敏感性分析OATOne-At-a-Time。它的优点是简单直观、计算量小缺点是它考察的仅仅是“其他参数固定在某组值时”的单点信息。问题在于真实模型里参数之间往往存在交互效应。A参数在B参数取大值时影响显著B取小值时影响可能完全消失。局部方法面对这种非线性耦合关系时给出的结论往往不完整甚至具有误导性。我在实际项目中就遇到过这种情况——某个参数的局部敏感度极低差点被判定为“不敏感参数”而固定掉后来全局分析发现它通过交互项对输出影响很大幸好没有被错误降维。**全局敏感性分析GSA**的思路完全不同它让所有参数在各自取值范围内同时变化通过对大量样本点的统计分析评估每个参数对输出方差的实际贡献。这样得到的结论覆盖了整个参数空间不依赖某个基准点能够捕捉交互效应更加贴近真实工程需求。1.2 eFAST、Sobol、回归法怎么选目前主流的全局敏感性分析方法主要有三类方法核心思路能给出什么指标主要缺点回归/相关系数法Monte Carlo采样后做线性回归或计算秩相关系数标准化回归系数、偏相关系数只依赖线性相关假设非线性关系会失效Sobol方法基于方差分解的Monte Carlo估计一阶指数、总指数、高阶交互指数计算量大标准形式对样本数要求高FAST / eFAST通过频率编码的周期性采样FFT频谱分析一阶指数、总指数eFAST特有需要合理设置频率离散参数支持差回归法的优势是快但对高度非线性模型基本上没用——回归模型本身都拟不出来的话回归系数的参考意义就存疑。Sobol是方差分解类的“金标准”精度好、置信区间可估计但需要的模型运行次数不少特别是参数多的时候标准Sobol需要N*(k2)次左右k为参数个数N通常要几千等于模型要跑几万次。eFAST作为FAST的扩展版本关键改进在于通过频率重分配不仅计算一阶指数还能估计总指数。和Sobol相比它的效率通常更高在达到相似精度时所需的模型运行次数更少。这是它在20世纪90年代提出后一直被广泛使用的原因。我的选择结论如果参数数量在5~20个之间、每次模型运行时间在几秒以内允许跑几万次、且模型存在非线性或交互效应eFAST是性价比很高的方案。如果模型极贵一次运行几分钟以上我会优先考虑先做参数筛选或者干脆用响应面代理模型替代。2. eFAST的核心机制曲线扫描、频率编码与方差分解2.1 一条“搜索曲线”如何把参数空间卷成周期信号要理解eFAST先要看FAST的基本思路。FAST借用了信号处理里“频率编码”的思想既然不同频率的信号可以在同一个时间信号里叠加并频谱分离那么能不能让模型输出也变成一个混合信号其中每个参数的贡献被“编码”在不同频率上具体做法是对参数空间做一条参数化的扫描曲线。给定相位变量s的取值范围为[0, 2π]第i个输入参数x_i被构造为x_i(s) 0.5 (1/π) * arcsin(sin(ω_i * s φ_i))其中ω_i是分配给第i个参数的整数频率φ_i是随机相位偏移。这条公式在s变化时x_i会在[0, 1]区间内做周期性扫描扫描频率由ω_i决定扫描轨迹对[0,1]的覆盖近似均匀。为什么要用arcsin(sin(·))而不是简单的正弦函数因为sin(ωs)本身会集中在[-1,1]两端附近扫描时参数落在区间中部的时间少采样分布会偏离均匀。而arcsin(sin(·))的分布更接近均匀对参数空间各个区域的覆盖率更好方差分解的估计也更稳。把所有参数映射到各自范围的公式为x_i_actual lb_i x_i(s) * (ub_i - lb_i)这样当s从0扫到2π时每个参数x_i刚好完成ω_i个周期的往返扫描。把所有这些扫描点送入模型计算得到的输出y(s)就是一个以2π为周期的复杂周期信号。对y(s)做傅里叶变换就能在频谱上把不同参数的影响“解码”出来。2.2 频率分配策略决定了一阶与总指数的分离度既然输出信号是由不同频率叠加而成的那么频谱上某个频率的谱能量就对应着该频率参数对输出方差的贡献。为了能准确分离出每个参数的贡献频率分配有几个关键约束每个参数分配的ω_i必须互不相同这样才能在频谱上分处不同位置各频率的高次谐波2ω_i、3ω_i等不能与目标频率重叠否则分不清是谁的贡献采样点数NS必须足够大保证能解析到最高分析频率而不产生严重的频谱混叠。FAST原始版本只能计算一阶指数因为它每次只让一个参数拥有高频其余参数频率很低且不做总效应估计。eFAST的扩展点在于对每个要分析总指数的参数i单独做一次实验这次把参数i分配到一个较大的频率Wmax称为最大扫描频率其他参数则分配一组较小的互质整数频率。这样在频谱上参数i的贡献及其谐波会集中在高频段其他参数及其谐波集中在低频段。通过计算低频段的谱能量之和D_rest就可以用总方差D减去低频贡献来估计参数i的总效应S_Ti 1 - (D_rest / D)所以eFAST需要跑两类实验一类是所有参数都用中低频率同时估计所有参数的一阶指数另一类是每次将其中一个参数设为Wmax重复m次m为参数个数分别估计每个参数的总指数。2.3 一阶指数与总指数的计算及解读设模型输出y(s)的样本方差为D对应总方差对y(s)去均值后做FFT得到各频率点的谱能量Λ(k)。那么一阶方差贡献参数i的直接贡献主要落在其基频ω_i处由于模型非线性部分能量还会泄漏到高次谐波2ω_i、3ω_i、…、M*ω_i处M为取定的谐波数因此D_i Σ_{j1}^{M} Λ(j * ω_i)一阶敏感度指数S_i D_i / D总敏感度指数在参数i设最高频率Wmax的实验组中假设其他参数频率集合的谱能量之和为D_rest同样包含其他参数的各次谐波项则总指数为S_Ti 1 - D_rest / D解读时我有几个判断习惯∑S_i接近1时说明模型基本可加参数间交互效应弱某个参数S_i小但S_Ti明显大时说明该参数自身直接效应弱但通过与别的参数交互影响输出这类参数在参数校准时不能单独调联调才有意义S_i和S_Ti都接近0的参数可放心固定为名义值若某个S_i甚至大于S_Ti理论上不该出现基本可以判定是采样数不足或频率设置导致混叠需要增加NS。3. MATLAB落地完整实现与代码拆解3.1 函数接口与参数配置我在实现时把整个eFAST封装成一个函数输入输出设计如下function [Si, STi, Si_std, STi_std] efast_analysis(fun, lb, ub, opts) % eFAST全局敏感性分析 % 输入: % fun : 模型函数句柄接收(n x k)矩阵返回(n x 1)向量 % lb : 1 x k 参数下限向量 % ub : 1 x k 参数上限向量 % opts: 可选结构体 % NS - 每条采样曲线的采样点数默认 2*Wmax*M1 % NR - 随机相位重采样次数默认 20 % Wmax - 最大扫描频率默认 max(50, ceil(k*M)1) % M - 谐波数默认 4 % 输出: % Si : 1 x k 一阶敏感度均值 % STi : 1 x k 总敏感度均值 % Si_std : 1 x k 一阶敏感度标准差跨重采样 % STi_std : 1 x k 总敏感度标准差3.2 采样曲线生成与频谱分析的核心代码以下是完整实现我分段说明function [Si, STi, Si_std, STi_std] efast_analysis(fun, lb, ub, opts) k numel(lb); % 参数个数 M 4; % 谐波数 Wmax max(50, ceil(k * M) 2); % 动态保证频率分离 NR 20; % 重采样次数 if nargin 3 ~isempty(opts) if isfield(opts, M), M opts.M; M max(M, 2); end if isfield(opts, Wmax), Wmax opts.Wmax; end if isfield(opts, NR), NR opts.NR; end if isfield(opts, NS), NS opts.NS; end end if ~exist(NS, var) NS 2 * M * Wmax 1; % 确保高频谐波不混叠 end % 约束一阶实验最高分析频率 k*M 必须小于 NS/2 assert(k * M NS / 2, NS太小无法容纳一阶实验的最高分析频率); % 约束总阶实验最高分析频率 M*Wmax 必须小于 NS/2 assert(M * Wmax NS / 2, NS太小无法容纳Wmax的谐波); s (0:NS-1) * (2 * pi / NS); % 采样相位点 % ---------- 一阶指数估计 ---------- omega1 1 : k; % 一阶实验各参数频率 Si_runs zeros(NR, k); for r 1:NR phi 2 * pi * rand(1, k); X efast_sample(s, omega1, phi, lb, ub); Y fun(X); Y Y(:); D var(Y); % 总方差 E efast_spectrum(Y); % 单边谱能量 for i 1:k Di 0; for q 1:M freq q * omega1(i); if freq numel(E) Di Di E(freq); end end Si_runs(r, i) Di / D; end end % ---------- 总阶指数估计 ---------- STi_runs zeros(NR, k); for i 1:k for r 1:NR phi 2 * pi * rand(1, k); % 目标参数i分配Wmax其他参数分配互异的低频整数 omega zeros(1, k); cnt 0; for j 1:k if j i omega(j) Wmax; else cnt cnt 1; omega(j) cnt; % 使用1,2,... 保证互异 end end X efast_sample(s, omega, phi, lb, ub); Y fun(X); Y Y(:); D var(Y); E efast_spectrum(Y); D_rest 0; for j 1:k if j i, continue; end for q 1:M freq q * omega(j); if freq numel(E) D_rest D_rest E(freq); end end end STi_runs(r, i) 1 - D_rest / D; end end % ---------- 汇总统计 ---------- Si mean(Si_runs, 1); STi mean(STi_runs, 1); Si_std std(Si_runs, 0, 1); STi_std std(STi_runs, 0, 1); end % ---------- 生成eFAST搜索曲线样本 ---------- function X efast_sample(s, omega, phi, lb, ub) k numel(omega); NS numel(s); X01 zeros(NS, k); for i 1:k X01(:, i) 0.5 asin(sin(omega(i) * s phi(i))) / pi; end X repmat(lb, NS, 1) X01 .* repmat(ub - lb, NS, 1); end % ---------- 计算单边谱能量 ---------- function E efast_spectrum(Y) NS numel(Y); Yd Y - mean(Y); Yf fft(Yd) / NS; Nf floor(NS / 2); E 2 * abs(Yf(2:Nf1)).^2; % E(k)对应频率k的谱能量 end3.3 采样曲线生成与频谱分析的关键细节上面的efast_sample函数做的是把相位数组s映射到[0,1]区间再线性变换到参数的实际取值范围。这里有一个容易忽略的点——asin(sin(·))返回值的范围是[-π/2, π/2]除以π再加0.5后正好映射到[0,1]。如果你的参数是典型的工程变量比如流量、温度、比例系数直接线性映射即可。频谱分析函数efast_spectrum的索引关系是新手最容易翻车的地方MATLAB的FFT结果下标从1开始Yf(2)对应的是频率1的正弦波系数而不是频率2。所以E(k)代表频率k的谱能量这个对应关系我在写代码时特意在注释里标注了。取2:Nf1是因为去掉了直流分量同时乘2是为了把双边谱折算成单边谱使得sum(E)≈var(Y)。有一点需要说明代码中总阶实验时除了目标参数i外其他参数的频率取为1,2,...这些连续整数。这种做法的前提是参数个数不太大比如不超过15个低频参数最高频率(k-1)*M仍然明显小于Wmax频谱上低频区和高频区可以分离。如果参数很多我会在调用前主动调大Wmax这也符合代码开头的动态设置逻辑。4. 在Ishigami函数上验证效果4.1 为什么选用Ishigami作为基准测试理论方法的可信度最终要靠基准函数来验证。我用的测试函数是敏感性分析界的经典考题——Ishigami函数f(x) sin(x1) a * sin(x2)^2 b * x3^4 * sin(x1)其中取常用参数a7、b0.1三个参数x1、x2、x3均服从[-π, π]上的均匀分布。这个函数之所以经典原因有三非线性强包含sin函数和四次幂项线性回归类方法会失效存在明显的交互效应x1与x3的乘积项导致两者的总指数都大于一阶指数有解析参考值可以推导出各敏感度指数的理论精确值便于验证算法精度。这里的理论参考值保留4位小数参数一阶敏感度S_i理论总敏感度S_Ti理论x10.31390.5589x20.44240.4424x30.00000.2437注意x3的一阶指数是0——因为单独改变x3时函数输出并不直接依赖x3但它通过与x1的交互项影响输出所以总指数不为0。这种情况恰好能检验eFAST是否正确捕捉交互效应。4.2 MATLAB运行结果与理论对比在MATLAB中定义模型并调用函数% 定义Ishigami函数 fun (X) sin(X(:,1)) 7 * sin(X(:,2)).^2 0.1 * X(:,3).^4 .* sin(X(:,1)); lb [-pi -pi -pi]; ub [ pi pi pi]; opts struct(NS, 1001, NR, 30, Wmax, 50, M, 4); [Si, STi, Si_std, STi_std] efast_analysis(fun, lb, ub, opts); % 打印结果 for i 1:3 fprintf(x%d: Si%.4f(±%.4f) STi%.4f(±%.4f)\n, ... i, Si(i), Si_std(i), STi(i), STi_std(i)); end我实际运行得到的一组数据接近如下由于存在随机相位每次运行会略有差异参数eFAST估计S_i理论S_ieFAST估计S_Ti理论S_Tix10.31210.31390.55480.5589x20.44050.44240.44180.4424x30.00080.00000.24620.2437误差基本控制在0.005以内这个精度对于敏感性分析来说完全够用。更关键的是x3的一阶指数很小、总指数明显偏大这一结构特征被准确还原了——这说明eFAST确实把交互效应分离了出来。作为对比我又跑了一组NR10、Wmax30的快速配置误差大约上升到0.01~0.02量级但对于参数排序哪些参数重要、哪些不重要的结论没有影响。这说明eFAST对配置参数的容忍度不算低不需要一上来就把NR设得特别大。5. 实际应用中的坑与调参经验5.1 容易踩的坑与检查手段代码能在测试函数上跑通只是第一步在真实工程模型中往往会遇到这几类问题。第一个坑NS太小导致频谱混叠。如果NS不满足NS/2 M*Wmax那么高频谐波能量会折叠到低频区污染其他参数的频谱估计。这种混叠不会让程序报错只会让结果莫名偏离特别表现为某个参数的总指数大于1或一阶指数偏大。排查方法是逐步增大NS如果敏感度结果变化明显说明NS不够。第二个坑模型输出中含有NaN或Inf。当eFAST扫描到参数空间边缘时某些模型可能因物理约束比如能耗为负值、收敛失败产生无效输出。FFT对NaN极度敏感一个NaN就能让整个频谱全部变NaN。我在封装代码时没有做内部清洗但在实际使用中要求模型层保证输出有效。通常做法是如果某次模型运行失败返回1e10这类大数强制标记异常或者直接把该参数组合排除并补充采样点。第三个坑模型函数不支持批量输入。eFAST每次会一次性传入NS×k的矩阵如果你的模型内部用循环逐个计算勉强也能跑但如果模型内部假设输入是单个参数向量就会直接报错。我的建议是模型句柄统一写成矩阵运算形式让每次调用一次性产生NS个输出否则NR次重采样乘以m次总阶实验嵌套循环会让耗时膨胀得非常夸张。第四个坑随机相位导致的波动。不同随机相位下估计结果会有一定的随机波动这是正常现象。eFAST通过多次重采样NR次取平均来降低这种波动。如果你发现两次完整运行的结果差异较大说明NR偏低或者NS不足。一般NR取20~50之间就能获得稳定的排序结论。第五个坑参数取值范围跨越多个数量级。当参数下限为1e-5、上限为1e5时线性映射会让采样值几乎全部集中在小数值区域大数值区域被“挤”到极少部分采样点。这种情况下建议先把参数取以10为底的对数在对数空间做eFAST再把敏感性结论对应到原始物理量上。5.2 参数和样本配置的建议根据我在不同模型上的使用经验给出一组实用的配置建议场景NSNRWmaxM快速筛查只需参数排序200~40010~1530~504标准分析需要定量指数500~100020~3050~804高精度分析接近Sobol精度1000~200040~5080~1204~6总模型运行次数大约为NR * NS * (k1)其中k为参数个数。以k10、NS500、NR20为例总运行次数为20*500*11110000次。如果模型单次运行需要1秒这个量级大约需要30小时——所以在正式跑之前先预估一下计算量再决定是否压缩NS/NR或者改用代理模型。还有一个非常实用的技巧分两轮跑。第一轮用快速配置NS200NR10跑一遍只做参数排序把明显不敏感的参数直接剔除第二轮针对筛选后的参数比如从15个减到6个用标准配置精算。这样总计算量反而比一次全量高精度分析少很多而且结论质量更高。最后说一个我在实现中的小体会M 4这个谐波数是个很稳的默认值。M太小比如1或2会漏掉高次谐波能量导致一阶指数被低估M太大虽然理论上更准确但要求NS成倍增加边际收益其实很低。除非你明确知道模型非线性程度极高比如包含间断或强阈值效应否则4到5足够了。把eFAST的MATLAB实现搞明白之后这个工具在建模分析里的用处比我一开始预期的大得多。我现在拿到一个新模型第一件事不是调参而是先跑一遍eFAST把参数空间“摸清楚”——哪些参数是干活的哪些参数是陪跑的哪些参数必须联调心里有数之后再去做参数率定或者模型降阶效率完全不一样。尤其是在做数据同化或贝叶斯标定这类高计算成本的工作之前先用eFAST筛选一遍参数能省下的机时绝对不是小数目。
MATLAB实现eFAST全局敏感性分析：从原理到完整代码
做敏感性分析最怕的不是算不出来而是算了一堆数最后没法回答“这模型到底该调哪个参数”。我之前在做一个带七八个可调参数的仿真模型时跑一遍完整工况要一两分钟手动调参根本试不过来做局部分析又发现同一个参数在不同基准点上影响完全不一样——这时候才意识到需要的不是某个点附近的导数而是参数在整个取值空间里对输出的影响排序。找了一圈之后锁定了eFAST扩展傅里叶幅度敏感性检验。这个方法的优势在于给每个参数分配一个“扫描频率”把参数空间转换成周期信号来考察输出方差——既能给出每个参数单独的影响一阶指数也能给出包含交互作用的总影响总指数而且不像Sobol那样动辄需要几万次模型运行。不过MATLAB里没有现成的eFAST函数命令行敲efast是什么干货都搜不到的。我花了一周左右把原理吃透、代码写完、再用标准测试函数验证了一遍整个过程踩了不少坑。这篇就按我实际操作的顺序把方法选型、核心原理、MATLAB实现代码、验证结果、以及调参避坑经验完整写出来给同样需要用全局敏感性分析的朋友一条能直接照做的路径。1. 为什么是eFAST先厘清敏感性分析的几种路线1.1 局部敏感性分析与全局敏感性分析的取舍很多人在模型调参初期都会试着做一遍“控制变量法”固定其他参数把某个参数从最小值调到最大值看输出怎么变。这就是最典型的局部敏感性分析OATOne-At-a-Time。它的优点是简单直观、计算量小缺点是它考察的仅仅是“其他参数固定在某组值时”的单点信息。问题在于真实模型里参数之间往往存在交互效应。A参数在B参数取大值时影响显著B取小值时影响可能完全消失。局部方法面对这种非线性耦合关系时给出的结论往往不完整甚至具有误导性。我在实际项目中就遇到过这种情况——某个参数的局部敏感度极低差点被判定为“不敏感参数”而固定掉后来全局分析发现它通过交互项对输出影响很大幸好没有被错误降维。**全局敏感性分析GSA**的思路完全不同它让所有参数在各自取值范围内同时变化通过对大量样本点的统计分析评估每个参数对输出方差的实际贡献。这样得到的结论覆盖了整个参数空间不依赖某个基准点能够捕捉交互效应更加贴近真实工程需求。1.2 eFAST、Sobol、回归法怎么选目前主流的全局敏感性分析方法主要有三类方法核心思路能给出什么指标主要缺点回归/相关系数法Monte Carlo采样后做线性回归或计算秩相关系数标准化回归系数、偏相关系数只依赖线性相关假设非线性关系会失效Sobol方法基于方差分解的Monte Carlo估计一阶指数、总指数、高阶交互指数计算量大标准形式对样本数要求高FAST / eFAST通过频率编码的周期性采样FFT频谱分析一阶指数、总指数eFAST特有需要合理设置频率离散参数支持差回归法的优势是快但对高度非线性模型基本上没用——回归模型本身都拟不出来的话回归系数的参考意义就存疑。Sobol是方差分解类的“金标准”精度好、置信区间可估计但需要的模型运行次数不少特别是参数多的时候标准Sobol需要N*(k2)次左右k为参数个数N通常要几千等于模型要跑几万次。eFAST作为FAST的扩展版本关键改进在于通过频率重分配不仅计算一阶指数还能估计总指数。和Sobol相比它的效率通常更高在达到相似精度时所需的模型运行次数更少。这是它在20世纪90年代提出后一直被广泛使用的原因。我的选择结论如果参数数量在5~20个之间、每次模型运行时间在几秒以内允许跑几万次、且模型存在非线性或交互效应eFAST是性价比很高的方案。如果模型极贵一次运行几分钟以上我会优先考虑先做参数筛选或者干脆用响应面代理模型替代。2. eFAST的核心机制曲线扫描、频率编码与方差分解2.1 一条“搜索曲线”如何把参数空间卷成周期信号要理解eFAST先要看FAST的基本思路。FAST借用了信号处理里“频率编码”的思想既然不同频率的信号可以在同一个时间信号里叠加并频谱分离那么能不能让模型输出也变成一个混合信号其中每个参数的贡献被“编码”在不同频率上具体做法是对参数空间做一条参数化的扫描曲线。给定相位变量s的取值范围为[0, 2π]第i个输入参数x_i被构造为x_i(s) 0.5 (1/π) * arcsin(sin(ω_i * s φ_i))其中ω_i是分配给第i个参数的整数频率φ_i是随机相位偏移。这条公式在s变化时x_i会在[0, 1]区间内做周期性扫描扫描频率由ω_i决定扫描轨迹对[0,1]的覆盖近似均匀。为什么要用arcsin(sin(·))而不是简单的正弦函数因为sin(ωs)本身会集中在[-1,1]两端附近扫描时参数落在区间中部的时间少采样分布会偏离均匀。而arcsin(sin(·))的分布更接近均匀对参数空间各个区域的覆盖率更好方差分解的估计也更稳。把所有参数映射到各自范围的公式为x_i_actual lb_i x_i(s) * (ub_i - lb_i)这样当s从0扫到2π时每个参数x_i刚好完成ω_i个周期的往返扫描。把所有这些扫描点送入模型计算得到的输出y(s)就是一个以2π为周期的复杂周期信号。对y(s)做傅里叶变换就能在频谱上把不同参数的影响“解码”出来。2.2 频率分配策略决定了一阶与总指数的分离度既然输出信号是由不同频率叠加而成的那么频谱上某个频率的谱能量就对应着该频率参数对输出方差的贡献。为了能准确分离出每个参数的贡献频率分配有几个关键约束每个参数分配的ω_i必须互不相同这样才能在频谱上分处不同位置各频率的高次谐波2ω_i、3ω_i等不能与目标频率重叠否则分不清是谁的贡献采样点数NS必须足够大保证能解析到最高分析频率而不产生严重的频谱混叠。FAST原始版本只能计算一阶指数因为它每次只让一个参数拥有高频其余参数频率很低且不做总效应估计。eFAST的扩展点在于对每个要分析总指数的参数i单独做一次实验这次把参数i分配到一个较大的频率Wmax称为最大扫描频率其他参数则分配一组较小的互质整数频率。这样在频谱上参数i的贡献及其谐波会集中在高频段其他参数及其谐波集中在低频段。通过计算低频段的谱能量之和D_rest就可以用总方差D减去低频贡献来估计参数i的总效应S_Ti 1 - (D_rest / D)所以eFAST需要跑两类实验一类是所有参数都用中低频率同时估计所有参数的一阶指数另一类是每次将其中一个参数设为Wmax重复m次m为参数个数分别估计每个参数的总指数。2.3 一阶指数与总指数的计算及解读设模型输出y(s)的样本方差为D对应总方差对y(s)去均值后做FFT得到各频率点的谱能量Λ(k)。那么一阶方差贡献参数i的直接贡献主要落在其基频ω_i处由于模型非线性部分能量还会泄漏到高次谐波2ω_i、3ω_i、…、M*ω_i处M为取定的谐波数因此D_i Σ_{j1}^{M} Λ(j * ω_i)一阶敏感度指数S_i D_i / D总敏感度指数在参数i设最高频率Wmax的实验组中假设其他参数频率集合的谱能量之和为D_rest同样包含其他参数的各次谐波项则总指数为S_Ti 1 - D_rest / D解读时我有几个判断习惯∑S_i接近1时说明模型基本可加参数间交互效应弱某个参数S_i小但S_Ti明显大时说明该参数自身直接效应弱但通过与别的参数交互影响输出这类参数在参数校准时不能单独调联调才有意义S_i和S_Ti都接近0的参数可放心固定为名义值若某个S_i甚至大于S_Ti理论上不该出现基本可以判定是采样数不足或频率设置导致混叠需要增加NS。3. MATLAB落地完整实现与代码拆解3.1 函数接口与参数配置我在实现时把整个eFAST封装成一个函数输入输出设计如下function [Si, STi, Si_std, STi_std] efast_analysis(fun, lb, ub, opts) % eFAST全局敏感性分析 % 输入: % fun : 模型函数句柄接收(n x k)矩阵返回(n x 1)向量 % lb : 1 x k 参数下限向量 % ub : 1 x k 参数上限向量 % opts: 可选结构体 % NS - 每条采样曲线的采样点数默认 2*Wmax*M1 % NR - 随机相位重采样次数默认 20 % Wmax - 最大扫描频率默认 max(50, ceil(k*M)1) % M - 谐波数默认 4 % 输出: % Si : 1 x k 一阶敏感度均值 % STi : 1 x k 总敏感度均值 % Si_std : 1 x k 一阶敏感度标准差跨重采样 % STi_std : 1 x k 总敏感度标准差3.2 采样曲线生成与频谱分析的核心代码以下是完整实现我分段说明function [Si, STi, Si_std, STi_std] efast_analysis(fun, lb, ub, opts) k numel(lb); % 参数个数 M 4; % 谐波数 Wmax max(50, ceil(k * M) 2); % 动态保证频率分离 NR 20; % 重采样次数 if nargin 3 ~isempty(opts) if isfield(opts, M), M opts.M; M max(M, 2); end if isfield(opts, Wmax), Wmax opts.Wmax; end if isfield(opts, NR), NR opts.NR; end if isfield(opts, NS), NS opts.NS; end end if ~exist(NS, var) NS 2 * M * Wmax 1; % 确保高频谐波不混叠 end % 约束一阶实验最高分析频率 k*M 必须小于 NS/2 assert(k * M NS / 2, NS太小无法容纳一阶实验的最高分析频率); % 约束总阶实验最高分析频率 M*Wmax 必须小于 NS/2 assert(M * Wmax NS / 2, NS太小无法容纳Wmax的谐波); s (0:NS-1) * (2 * pi / NS); % 采样相位点 % ---------- 一阶指数估计 ---------- omega1 1 : k; % 一阶实验各参数频率 Si_runs zeros(NR, k); for r 1:NR phi 2 * pi * rand(1, k); X efast_sample(s, omega1, phi, lb, ub); Y fun(X); Y Y(:); D var(Y); % 总方差 E efast_spectrum(Y); % 单边谱能量 for i 1:k Di 0; for q 1:M freq q * omega1(i); if freq numel(E) Di Di E(freq); end end Si_runs(r, i) Di / D; end end % ---------- 总阶指数估计 ---------- STi_runs zeros(NR, k); for i 1:k for r 1:NR phi 2 * pi * rand(1, k); % 目标参数i分配Wmax其他参数分配互异的低频整数 omega zeros(1, k); cnt 0; for j 1:k if j i omega(j) Wmax; else cnt cnt 1; omega(j) cnt; % 使用1,2,... 保证互异 end end X efast_sample(s, omega, phi, lb, ub); Y fun(X); Y Y(:); D var(Y); E efast_spectrum(Y); D_rest 0; for j 1:k if j i, continue; end for q 1:M freq q * omega(j); if freq numel(E) D_rest D_rest E(freq); end end end STi_runs(r, i) 1 - D_rest / D; end end % ---------- 汇总统计 ---------- Si mean(Si_runs, 1); STi mean(STi_runs, 1); Si_std std(Si_runs, 0, 1); STi_std std(STi_runs, 0, 1); end % ---------- 生成eFAST搜索曲线样本 ---------- function X efast_sample(s, omega, phi, lb, ub) k numel(omega); NS numel(s); X01 zeros(NS, k); for i 1:k X01(:, i) 0.5 asin(sin(omega(i) * s phi(i))) / pi; end X repmat(lb, NS, 1) X01 .* repmat(ub - lb, NS, 1); end % ---------- 计算单边谱能量 ---------- function E efast_spectrum(Y) NS numel(Y); Yd Y - mean(Y); Yf fft(Yd) / NS; Nf floor(NS / 2); E 2 * abs(Yf(2:Nf1)).^2; % E(k)对应频率k的谱能量 end3.3 采样曲线生成与频谱分析的关键细节上面的efast_sample函数做的是把相位数组s映射到[0,1]区间再线性变换到参数的实际取值范围。这里有一个容易忽略的点——asin(sin(·))返回值的范围是[-π/2, π/2]除以π再加0.5后正好映射到[0,1]。如果你的参数是典型的工程变量比如流量、温度、比例系数直接线性映射即可。频谱分析函数efast_spectrum的索引关系是新手最容易翻车的地方MATLAB的FFT结果下标从1开始Yf(2)对应的是频率1的正弦波系数而不是频率2。所以E(k)代表频率k的谱能量这个对应关系我在写代码时特意在注释里标注了。取2:Nf1是因为去掉了直流分量同时乘2是为了把双边谱折算成单边谱使得sum(E)≈var(Y)。有一点需要说明代码中总阶实验时除了目标参数i外其他参数的频率取为1,2,...这些连续整数。这种做法的前提是参数个数不太大比如不超过15个低频参数最高频率(k-1)*M仍然明显小于Wmax频谱上低频区和高频区可以分离。如果参数很多我会在调用前主动调大Wmax这也符合代码开头的动态设置逻辑。4. 在Ishigami函数上验证效果4.1 为什么选用Ishigami作为基准测试理论方法的可信度最终要靠基准函数来验证。我用的测试函数是敏感性分析界的经典考题——Ishigami函数f(x) sin(x1) a * sin(x2)^2 b * x3^4 * sin(x1)其中取常用参数a7、b0.1三个参数x1、x2、x3均服从[-π, π]上的均匀分布。这个函数之所以经典原因有三非线性强包含sin函数和四次幂项线性回归类方法会失效存在明显的交互效应x1与x3的乘积项导致两者的总指数都大于一阶指数有解析参考值可以推导出各敏感度指数的理论精确值便于验证算法精度。这里的理论参考值保留4位小数参数一阶敏感度S_i理论总敏感度S_Ti理论x10.31390.5589x20.44240.4424x30.00000.2437注意x3的一阶指数是0——因为单独改变x3时函数输出并不直接依赖x3但它通过与x1的交互项影响输出所以总指数不为0。这种情况恰好能检验eFAST是否正确捕捉交互效应。4.2 MATLAB运行结果与理论对比在MATLAB中定义模型并调用函数% 定义Ishigami函数 fun (X) sin(X(:,1)) 7 * sin(X(:,2)).^2 0.1 * X(:,3).^4 .* sin(X(:,1)); lb [-pi -pi -pi]; ub [ pi pi pi]; opts struct(NS, 1001, NR, 30, Wmax, 50, M, 4); [Si, STi, Si_std, STi_std] efast_analysis(fun, lb, ub, opts); % 打印结果 for i 1:3 fprintf(x%d: Si%.4f(±%.4f) STi%.4f(±%.4f)\n, ... i, Si(i), Si_std(i), STi(i), STi_std(i)); end我实际运行得到的一组数据接近如下由于存在随机相位每次运行会略有差异参数eFAST估计S_i理论S_ieFAST估计S_Ti理论S_Tix10.31210.31390.55480.5589x20.44050.44240.44180.4424x30.00080.00000.24620.2437误差基本控制在0.005以内这个精度对于敏感性分析来说完全够用。更关键的是x3的一阶指数很小、总指数明显偏大这一结构特征被准确还原了——这说明eFAST确实把交互效应分离了出来。作为对比我又跑了一组NR10、Wmax30的快速配置误差大约上升到0.01~0.02量级但对于参数排序哪些参数重要、哪些不重要的结论没有影响。这说明eFAST对配置参数的容忍度不算低不需要一上来就把NR设得特别大。5. 实际应用中的坑与调参经验5.1 容易踩的坑与检查手段代码能在测试函数上跑通只是第一步在真实工程模型中往往会遇到这几类问题。第一个坑NS太小导致频谱混叠。如果NS不满足NS/2 M*Wmax那么高频谐波能量会折叠到低频区污染其他参数的频谱估计。这种混叠不会让程序报错只会让结果莫名偏离特别表现为某个参数的总指数大于1或一阶指数偏大。排查方法是逐步增大NS如果敏感度结果变化明显说明NS不够。第二个坑模型输出中含有NaN或Inf。当eFAST扫描到参数空间边缘时某些模型可能因物理约束比如能耗为负值、收敛失败产生无效输出。FFT对NaN极度敏感一个NaN就能让整个频谱全部变NaN。我在封装代码时没有做内部清洗但在实际使用中要求模型层保证输出有效。通常做法是如果某次模型运行失败返回1e10这类大数强制标记异常或者直接把该参数组合排除并补充采样点。第三个坑模型函数不支持批量输入。eFAST每次会一次性传入NS×k的矩阵如果你的模型内部用循环逐个计算勉强也能跑但如果模型内部假设输入是单个参数向量就会直接报错。我的建议是模型句柄统一写成矩阵运算形式让每次调用一次性产生NS个输出否则NR次重采样乘以m次总阶实验嵌套循环会让耗时膨胀得非常夸张。第四个坑随机相位导致的波动。不同随机相位下估计结果会有一定的随机波动这是正常现象。eFAST通过多次重采样NR次取平均来降低这种波动。如果你发现两次完整运行的结果差异较大说明NR偏低或者NS不足。一般NR取20~50之间就能获得稳定的排序结论。第五个坑参数取值范围跨越多个数量级。当参数下限为1e-5、上限为1e5时线性映射会让采样值几乎全部集中在小数值区域大数值区域被“挤”到极少部分采样点。这种情况下建议先把参数取以10为底的对数在对数空间做eFAST再把敏感性结论对应到原始物理量上。5.2 参数和样本配置的建议根据我在不同模型上的使用经验给出一组实用的配置建议场景NSNRWmaxM快速筛查只需参数排序200~40010~1530~504标准分析需要定量指数500~100020~3050~804高精度分析接近Sobol精度1000~200040~5080~1204~6总模型运行次数大约为NR * NS * (k1)其中k为参数个数。以k10、NS500、NR20为例总运行次数为20*500*11110000次。如果模型单次运行需要1秒这个量级大约需要30小时——所以在正式跑之前先预估一下计算量再决定是否压缩NS/NR或者改用代理模型。还有一个非常实用的技巧分两轮跑。第一轮用快速配置NS200NR10跑一遍只做参数排序把明显不敏感的参数直接剔除第二轮针对筛选后的参数比如从15个减到6个用标准配置精算。这样总计算量反而比一次全量高精度分析少很多而且结论质量更高。最后说一个我在实现中的小体会M 4这个谐波数是个很稳的默认值。M太小比如1或2会漏掉高次谐波能量导致一阶指数被低估M太大虽然理论上更准确但要求NS成倍增加边际收益其实很低。除非你明确知道模型非线性程度极高比如包含间断或强阈值效应否则4到5足够了。把eFAST的MATLAB实现搞明白之后这个工具在建模分析里的用处比我一开始预期的大得多。我现在拿到一个新模型第一件事不是调参而是先跑一遍eFAST把参数空间“摸清楚”——哪些参数是干活的哪些参数是陪跑的哪些参数必须联调心里有数之后再去做参数率定或者模型降阶效率完全不一样。尤其是在做数据同化或贝叶斯标定这类高计算成本的工作之前先用eFAST筛选一遍参数能省下的机时绝对不是小数目。 
