用广义分层抽样和随机仿真器加速结构性能化风险优化

性能化风险优化在结构工程里的核心矛盾一直很明确风险评估需要大量非线性时程分析优化又要在风险评估外层再套一层搜索两种计算成本叠在一起普通工作站很难在可接受时间内完成。针对这个问题“Stochastic Emulation using Generalized Stratified Sampling for Performance-Based Risk Optimization of Structures”这条技术路线提供了一个通用解法先用广义分层抽样Generalized Stratified Sampling简称 GSS挑出少量高信息量的结构分析样本再用随机仿真器Stochastic Emulator代替真实结构分析最后把风险指标放进优化循环寻找最优设计。本文不讨论论文原文的具体算例而是按这条技术主线拆解概念、给出一个最小 Python 示例、说明参数取舍、验证方式和排查路径。适合结构工程方向的研究生、做不确定性量化和概率设计的工程师以及想了解代理模型如何落地到工程优化中的开发者。1. 先理解随机仿真、仿真器和风险优化是怎么串起来的1.1 性能化风险优化为什么计算量会失控性能化抗震设计的思路是把结构性能目标从“不倒塌”扩展为“不同地震水平下满足不同位移或损伤指标”。放到风险评估里通常要计算某个设计方案在未来一定年限内的年超越概率或者期望年损失。计算流程大致是选一组设计参数例如基本周期、屈服强度、构件截面尺寸。在多个地震动强度水平下输入多条地面运动记录。每一条记录做一次非线性时程分析得到最大层间位移角、残余位移等响应。把响应与损伤状态和损失函数挂钩积分得到风险指标。问题在于一次时程分析可能是分钟到小时级一个设计点要跑几百条记录才能估计它的超越概率而优化算法又可能需要评估上千个设计点。如果直接套两层循环总计算量会变成“设计点数量 × 地震动数量 × 单次分析耗时”这正是性能化风险优化难以工程化的根本原因。1.2 Emulator 解决的是“评估次数”问题Emulator 在本文语境中指的是“随机仿真器”是一种用统计学习模型构造的代理函数。它的输入是结构设计参数和随机变量输出是结构响应或风险指标的预测值并且预测时不仅要给出均值还要给出不确定性范围。工程中常见的叫法是“代理模型surrogate model”但 emulator 更强调两点第一它学习的对象是随机仿真输出本身带有噪声第二它必须提供预测方差这样后续优化才能判断某个设计点到底是“确定地差”还是“因为训练样本少而看不清楚”。高斯过程回归Gaussian Process Regression天然满足这两个要求因此成为最常用的选择。相比直接在每个设计点做时程分析训练好的 GP 每次预测只需毫秒级时间可以让优化算法在同样的时间预算里多评估几个数量级的候选方案。1.3 GSS 解决的是“训练样本”问题仿真器再快训练数据仍然来自真实结构分析所以样本怎么采直接决定训练成本和精度。广义分层抽样要回答的问题就是在有限分析预算下样本点放在参数空间的哪些位置才能让仿真器对风险目标的估计误差最小。经典分层抽样的做法是把输入空间按概率分成若干层每层按比例分配样本。广义分层抽样则把“按比例”扩展成“按风险贡献分配”层边界可以根据物理意义或者变量分位数确定层内样本量可以根据该层对目标响应方差的贡献来调整也就是类似 Neyman 分配的思路。层内还可以再用拉丁超立方或者简单随机抽样填充这样既保证样本覆盖整个参数空间又不会把预算浪费在对风险影响很小的区域。当一个结构设计变量的微小变化会导致响应剧烈变化时普通随机抽样很容易漏掉这些敏感区域而分层结构能显著降低风险估计的方差。1.4 完整链路从样本到优化再到校核把上面的概念串起来整个方法的主线是定义设计变量、随机变量和风险指标。用 GSS 生成一组覆盖参数空间的分析样本。对每个样本执行真实结构分析得到响应或风险输出。用这些输入输出对训练随机仿真器。在优化循环中用仿真器替代真实分析快速评估目标函数。对优化结果用少量真实分析校核。这条链路上GSS 决定数据质量仿真器决定计算效率优化器决定搜索效率。任一个环节偏弱最终设计都不可信所以后面几节分别展开。2. 方法框架从设计变量到风险目标2.1 输入、随机变量和风险输出要一开始就分清楚很多失败的实践问题不是算法本身而是变量分类不清。至少要区分三类对象设计变量优化算法可以改变的量例如基本周期、强度折减系数、构件截面尺寸。随机变量无法直接控制或者只能控制分布的量例如地震动强度、地震动记录之间的变异、材料强度离散性。风险输出用于评估方案优劣的量例如最大层间位移角的超越概率、倒塌概率、期望年损失。这三类对象一旦混在一起后面的采样和优化都会失去意义。设计变量要进优化器的边界随机变量要进采样器的分层范围风险输出是目标函数。训练仿真器时输入列应同时包含设计变量和随机变量因为我们的目标是用它们共同预测输出。2.2 广义分层抽样的设计步骤GSS 的落地可以按四步进行。第一步确定分层维度。一般选择对响应影响最大的 2 到 5 个变量做分层变量太多会导致单元数量爆炸。比如地震动强度、结构周期、强度折减系数就是常见的高影响维度。第二步确定层边界。可以用等间距、分位数或者根据物理阈值确定。分位数方式能保证每个边界内的样本出现概率大致一致适合重尾分布等间距方式实现简单适合均匀设计空间。第三步确定每层样本量。最简单是体积比例分配更高效的做法是用 Neyman 分配即样本量与层内标准差成正比n_h n_total * (w_h * sigma_h) / sum(w_i * sigma_i)其中 w_h 是第 h 层的权重sigma_h 是第 h 层响应的标准差。问题是 sigma_h 在分析前未知所以工程上常用两阶段策略先抽少量样本估计每层方差再分配剩余预算。第四步在每一层内部用拉丁超立方或者均匀随机采样填充。层内填充方式决定了样本的空间散布质量拉丁超立方通常优于纯随机。2.3 仿真器训练与噪声建模随机仿真的输出有噪声。即使用同一个设计参数换一条地震动记录响应也不一样。因此 GP 的核函数里必须包含噪声项否则模型会误把噪声当成信号预测方差偏小后续优化会过度相信某些假象。训练时要把数据划分成训练集和测试集用独立的测试样本评估 R²、RMSE 和 2σ 覆盖比例。覆盖比例的含义是真实响应落在预测均值 ±2 倍标准差范围内的样本占比。理想情况下这个比例应该接近 95%如果明显偏低说明 GP 的不确定性程度没有校准好不能立即用于优化。2.4 优化循环的设计优化目标可以直接用风险指标本身也可以把风险指标作为约束把建造成本作为目标。典型写法是目标函数最小化年超越概率或者最小化期望年损失。约束条件位移角响应不超限、构件强度不失效、可施工性边界。由于 GP 能输出预测标准差优化时可以区分“确定的好”和“不确定的好”。工程上通常取预测分布的某个高分位点作为保守评估值例如用均值加 1.5 倍标准差。这样即使仿真器在某个区域有误差优化结果也不会因为噪声而被带偏。3. 最小可复现示例Python GP 仿真器 简化结构黑盒3.1 环境准备示例代码依赖以下 Python 包pip install numpy scipy scikit-learn matplotlib版本方面没有特别严格的限制使用近两年发布的 numpy、scipy 和 scikit-learn 版本即可。示例中用于演示的“真实结构分析”是一个简化单自由度滞回模型实际项目里应当替换为 OpenSees、ABAQUS 或 PERFORM-3D 的分析调用但接口保持一致后面的采样、训练、优化代码可以复用。3.2 定义一个黑盒结构响应函数下面代码用 Bouc-Wen 平滑滞回模型模拟一个单自由度结构在地震动作用下的最大位移角。这里的关键不是模型本身而是它扮演“昂贵真实分析”的角色。import numpy as np from scipy.integrate import odeint def synthetic_ground_motion(pga, duration20.0, dt0.005, seed7): 生成一条人工地震动时程仅用于方法演示。 rng np.random.default_rng(seed) t np.arange(0.0, duration, dt) noise rng.standard_normal(len(t)) smooth np.convolve(noise, np.ones(40) / 40, modesame) envelope np.exp(-2.5 * t / duration) acc smooth * envelope acc acc / (np.max(np.abs(acc)) 1e-12) * pga return t, acc def run_sdof(theta, im_record): Bouc-Wen 单自由度滞回模型返回按屈服位移归一化的最大位移。 theta [T, ksi, ry] T : 结构基本周期 ksi : 阻尼比 ry : 屈服强度折减系数 实际工程中这个函数应替换为 OpenSees / ABAQUS 的分析调用。 T, ksi, ry theta t, ag im_record omega 2.0 * np.pi / T m 1.0 k m * omega ** 2 c 2.0 * m * omega * ksi F_y (0.2 * 9.81) / ry u_y F_y / k alpha 0.05 A, beta, gamma, n 1.0, 0.5, 0.5, 1.0 def deriv(s, tt): u, v, z s a_g np.interp(tt, t, ag) F alpha * k * u (1.0 - alpha) * k * u_y * z u_dd -a_g - (c * v F) / m z_d (A * v - beta * np.abs(v) * np.abs(z) ** (n - 1) * z - gamma * v * np.abs(z) ** n) / u_y return [v, u_dd, z_d] sol odeint(deriv, [0.0, 0.0, 0.0], t) u sol[:, 0] return float(np.max(np.abs(u)) / u_y)这段代码只用于说明接口设计参数不代表真实结构正式项目中的分析模型必须经过验证。理解重点在于run_sdof的输入是一个设计参数数组和一个地震动输出是一个标量响应真实有限元模型也可以用同样的函数封装。3.3 用广义分层抽样生成训练样本下面的采样函数把每个维度均匀切成若干层然后按体积占比分配样本数层内做均匀随机采样。这里实现的是最简单的体积比例分层实际项目中可以在这个基础上加入方差权重。def gss_sample(bounds, n_total, n_strata4, seed42): 在参数空间上做分层采样。 bounds : 每个变量的 [下界, 上界] 列表 n_total : 期望样本总量 n_strata: 每个维度切分的层数 rng np.random.default_rng(seed) dim len(bounds) widths [hi - lo for lo, hi in bounds] edges [np.linspace(lo, hi, n_strata 1) for lo, hi in bounds] points [] for idx in np.ndindex(*([n_strata] * dim)): vol_frac 1.0 for d in range(dim): vol_frac * (edges[d][idx[d] 1] - edges[d][idx[d]]) / widths[d] n_cell max(1, int(round(n_total * vol_frac))) for _ in range(n_cell): p [rng.uniform(edges[d][idx[d]], edges[d][idx[d] 1]) for d in range(dim)] points.append(p) return np.array(points)使用方式bounds [(0.5, 3.0), (0.02, 0.08), (2.0, 8.0)] # T, ksi, ry X gss_sample(bounds, n_total80, n_strata4, seed42) y [] for theta in X: record_seed int(abs(np.sum(theta)) * 100) % 10000 im synthetic_ground_motion(pga0.4, seedrecord_seed) y.append(run_sdof(theta, im)) y np.array(y)这里每个样本使用不同 seed 生成地震动模拟记录间变异性。样本量 80 是为了让示例快速运行正式项目中样本量可能要上百甚至更多。3.4 训练高斯过程仿真器GP 的核函数采用常数核乘以 RBF 核再加白噪声核。白噪声核对应随机仿真的噪声项不能省略。from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel from sklearn.model_selection import train_test_split from sklearn.metrics import r2_score, mean_squared_error X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.25, random_state1 ) kernel (ConstantKernel(1.0) * RBF(length_scale[0.3, 0.3, 0.3]) WhiteKernel(1e-3)) gp GaussianProcessRegressor( kernelkernel, n_restarts_optimizer5, normalize_yTrue ) gp.fit(X_train, y_train) y_pred, y_std gp.predict(X_test, return_stdTrue) print(R2 , r2_score(y_test, y_pred)) print(RMSE , mean_squared_error(y_test, y_pred, squaredFalse)) within np.mean(np.abs(y_test - y_pred) 2.0 * y_std) print(2σ 覆盖比例 , within)normalize_yTrue让模型先标准化输出再拟合对响应量级差异大的问题很有帮助。n_restarts_optimizer表示从多个起始点优化超参数避免陷入局部最优。3.5 风险优化循环优化目标使用一个简化风险指标GP 预测分布超过位移角阈值的概率。这个指标不是真正的年超越概率但足以演示“用带不确定性的预测做优化”的写法。from scipy.optimize import differential_evolution from math import erf P_THRESHOLD 0.02 # 位移角阈值示例用 def risk_objective(theta): theta np.asarray(theta).reshape(1, -1) mean, std gp.predict(theta, return_stdTrue) z (P_THRESHOLD - mean[0]) / (std[0] * np.sqrt(2.0)) return 0.5 * (1.0 - erf(z)) res differential_evolution( risk_objective, bounds, seed42, tol1e-8, polishTrue ) print(最优设计参数 , res.x) print(风险指标 , res.fun)differential_evolution是 scipy 内置的差分进化算法适合维度不高、边界清楚的优化问题。polishTrue会在差分进化结束后用局部优化器进一步精修解。3.6 运行与预期输出按顺序运行完整代码预期可以看到R² 在 0.85 以上说明仿真器对测试集预测能力可用。RMSE 小于响应标准差的一半说明均方误差控制得较好。2σ 覆盖比例在 0.90 到 0.98 之间说明不确定性范围基本合理。优化输出一组满足边界条件的设计参数风险指标值明显低于随机初始设计。如果 R² 低于 0.7不要急着调优化器先检查采样覆盖、样本量和核函数设置。4. 关键参数怎么选先看含义再看影响4.1 样本量、分层数与层边界样本量决定训练成本。GP 在小样本下表现好但样本太少会漏掉参数空间的敏感区域。一般先跑 50 到 100 个样本如果测试集精度不足再增加。不要一开始就追求上千样本因为真实结构分析成本很高。分层数 n_strata 取决于每个维度的非线性程度。对响应平滑的参数3 层足够响应在某段会剧烈变化的参数可以增加到 5 层。每维超过 5 层在小样本下反而会导致单个单元内样本太少方差估计不稳定。层边界的确定要结合物理含义。例如地震动强度按分位数切分可以保证每个强度区间都有代表样本结构周期按等间距切分则便于后续解释。不要盲目用均匀网格覆盖全空间因为很多参数空间本身存在明显的高概率区间和低概率区间。4.2 仿真器结构与核函数对大多数结构响应问题RBF 核加白噪声核就够用。RBF 核的 length_scale 控制函数平滑程度可以由优化器自动估计不必手动调整。如果响应存在明显的线性趋势可以加一个线性核如果响应维度很高则要考虑降维或者改用稀疏高斯过程。随机仿真的重复性噪声量级很重要。WhiteKernel 的初始值设置过小模型会把噪声当信号设置过大模型会过度平滑。建议先用几组完全相同的输入重复分析直接估计响应噪声量级再设置 WhiteKernel 的初始值。4.3 优化器与风险指标优化器选择取决于目标函数形态。差分进化适合非光滑、多峰问题但存在边界附近的样本可能不稳定贝叶斯优化适合分析成本更高的场景它利用 GP 的预测方差决定下一步在哪里采样。示例代码使用差分进化工程中往往先做一次差分进化搜索再用局部优化器精修。风险指标的定义会直接影响最优设计。如果只用响应均值优化结果可能低估风险如果采用均值加 1.5 到 2 倍标准差则结果偏保守但更稳。建议在优化前至少试两种风险代理对比最优设计差异差异过大说明仿真器不确定性还没有收敛。关键参数速查表如下参数含义常见取值调大的影响调小的影响n_total训练样本总量50 ~ 200精度更高分析成本上升可能漏掉敏感区n_strata每维分层数3 ~ 5覆盖更细单元样本变少覆盖粗糙RBF length_scaleGP 核长度尺度自动估计预测更平滑预测更抖动WhiteKernel noise噪声项由重复试验估计过度平滑过拟合噪声n_restarts_optimizer超参数优化重启次数3 ~ 10更稳定更耗时可能陷入局部最优风险分位系数优化用保守程度1.5 ~ 2.0结果更保守结果更激进5. 结果验证不能只跑通还要校核模型可信度5.1 仿真器精度验证指标训练集上的 R² 没有参考意义必须用独立测试集。三个指标一起看R² 衡量整体拟合比例0.85 以上算基本可用。RMSE 衡量绝对误差要结合响应量级判断例如最大位移角均值 0.03RMSE 0.003 就是可接受的。2σ 覆盖比例衡量不确定性校准理想约 0.95。低于 0.90 说明预测标准差偏小高于 0.99 说明预测标准差偏大。更严格的做法是 k 折交叉验证尤其是样本量小时。每折训练一次 GP统计所有测试点的误差分布。5.2 风险估计的方差对比仿真器精度不等于风险估计精度。要验证 GSS 是否真的比简单随机抽样好可以用相同训练样本量分别用简单随机抽样和 GSS 各抽一批样本训练两个仿真器再对同一个设计点估计超越概率。重复 20 到 50 次比较两种采样方式下超越概率估计的方差。GSS 的方差应当更小或者达到相同方差时所需样本更少。工程中可以把“样本量-风险估计方差”曲线画出来横轴是分析次数纵轴是风险指标的标准差。曲线越早进入低方差平台说明采样策略越高效。这个验证往往揭示一个重要事实提升风险估计稳定性主要靠采样策略而不只是靠增大 GP 的复杂度。5.3 优化结果必须用真实分析校核仿真是近似优化结果不能直接放行。对优化得到的最优设计点用真实结构分析重新评估风险指标如果真实值与仿真器预测偏差超过工程允许范围需要在该点附近补采样本重新训练仿真器再重新优化直到偏差收敛。这本质上是主动学习的思想。报告输出时至少包含样本设计表、训练集与测试集误差指标、风险估计方差对比、优化历史曲线、最优设计点真实分析与仿真器预测的对比。这样别人才能判断结果可信度。6. 常见问题与排查路径6.1 问题现象与处理方案问题现象常见原因检查方式处理建议GP 测试集 R² 很低训练样本覆盖不足或样本量偏小检查样本是否集中在局部区域增加 GSS 样本量或者提高该维度分层数预测标准差整体过大噪声项初始值偏大或数据噪声大查看 WhiteKernel 优化后的值用重复试验估计真实噪声修正核函数初始值2σ 覆盖比例远低于 0.90模型低估不确定性检查是否删除了噪声项恢复 WhiteKernel 或增大噪声初始值优化结果落在边界附近震荡风险指标在该区域不敏感打印优化历史增加惩罚项或改用局部-全局两阶段优化真实分析成本太高训练集生成太慢每个样本都跑大量地震动检查样本量是否过度先减少记录数量做预筛再用 GSS 集中预算分层单元内样本只有 1 到 2 个n_strata 过高而 n_total 不足查看各单元样本数降低每维分层数或者用拉丁超立方替代纯随机6.2 排查顺序遇到问题先按下面的顺序定位不要直接换算法检查输入变量是否归一化。GP 对尺度敏感量级差异过大的变量会导致核函数长度尺度估计失效。检查训练数据是否存在重复或异常值。重复样本会放大噪声项异常值会拉偏 RBF 核。检查采样覆盖。直接打印样本在每个维度的直方图看是否出现大段空白。