做结构抗震评估这行当迟早要碰增量动力分析IDA和易损性曲线。不管你是写论文、做韧性评估还是给既有建筑做地震风险研判都绕不开“结构在多大的地震动下会坏到什么程度”这个问题。我最早接触IDA是在做某框架结构的倒塌概率评估时被“打点、拟合、画曲线”这几个字折磨得不轻后来把Matlab代码一步步拆开才真正明白这套方法并不神秘难点全藏在细节里。这篇文章不打算讲教科书式的理论推导而是从实战角度把IDA到易损性曲线的完整流程、Matlab实现思路和我在实际计算中踩过的坑都摊开说清楚。适合正在写结构抗震方向论文的研究生、需要用概率方法评估结构性能的工程师以及那些想在Matlab里亲手把易损性曲线跑出来却不知道从哪下手的同学。1. IDA到底在算什么IM、EDP、极限状态三者的关系1.1 一句话说清增量动力分析增量动力分析Incremental Dynamic AnalysisIDA本质上是一种“对结构反复做非线性时程分析”的试验方法。差别在于每一次分析的地震动强度是逐步增大的然后记录下结构在不同强度下的动力响应。具体操作流程是这样选取一批地震动记录每一条记录都按不同的调幅系数进行缩放从小到大形成一系列强度水准。对每一个强度水准把缩放过后的地震波输入到结构模型里做一次完整的非线性时程分析提取出该工况下的最大响应值比如最大层间位移角。最后以地震动强度指标IM为横坐标、结构需求参数EDP为纵坐标把同一条地震动记录在不同强度下的响应点连成一条曲线多条记录对应多条曲线就得到了经典的IDA曲线族。这里的IMIntensity Measure可以是峰值加速度PGA、峰值速度PGV也可以是结构第一周期处的谱加速度Sa(T1, 5%)。EDPEngineering Demand Parameter通常取最大层间位移角、顶点位移、塑性转角等。之所以反复强调“同一个结构、同一条地震动、逐步增大强度”是因为IDA曲线刻画的是结构从弹性到屈服、再到倒塌的完整演化过程这是普通单次时程分析给不了的信息。1.2 易损性曲线是在回答什么问题结构易损性曲线本质上是一个概率表达当地震动强度等于某个IM值时结构达到或超过某一极限状态的概率是多少。比如当PGA为0.3g时最大层间位移角超过1/50的概率是15%这就是易损性曲线上的一个点。为什么不能直接给一个确定性的临界强度因为地震动有极强的随机性两条峰值加速度相同的记录实际对结构造成的影响可能差出一倍。结构自身的材料强度、刚度也有离散性。所以工程上更倾向于回答“在某个强度下有百分之多少的可能性进入某种破坏状态”而不是“这个强度一定破坏”。易损性曲线就是这种概率化回答的可视化结果。IDA和易损性曲线之间的桥梁是这样搭起来的IDA曲线提供了每个IM水平下结构响应的分布而极限状态定义了“破坏”的门槛。把两者一对照就能算出每个IM水平下的失效频率再通过统计拟合得到一条光滑的、连续的概率曲线。整个过程如果用Matlab实现代码架构并不复杂难在数据组织和非线性时程分析的稳定计算。2. Matlab代码的整体框架从地震波导入到曲线输出的数据流2.1 模块划分与数据组织思路写IDA代码最忌讳一上来就写一个巨大的主程序所有操作揉在一起。我建议按数据流把代码拆成五个模块每个模块负责一件事输出特定的中间数据。这样调试起来非常舒服哪一步出问题一眼就能定位。模块功能主要输入主要输出地震动处理模块读取、截断、基线校正、统一时间步长原始加速度时程文件标准化后的加速度时程数组缩放与IM计算模块计算PGA或Sa(T1)生成缩放系数标准化地震动、结构周期T1缩放后的时程、IM值时程分析模块执行非线性动力积分提取EDP缩放的时程、结构模型参数最大层间位移角等EDP数据汇总模块按IM从小到大排列所有工况结果各工况的IM与EDPIDA曲线的原始数据表易损性拟合模块计算失效概率并拟合对数正态分布极限状态阈值、IDA数据表易损性曲线参数与绘图数据这里要特别说一下数据组织方式。通常会设计一个二维或三维矩阵比如EDP_all(i, j)表示第i条地震动、第j个强度水平下的最大层间位移角。同时再保存一个IM_all(i, j)记录对应的强度值。这样后面画IDA曲线、统计失效概率、拟合易损性曲线时直接取矩阵里的列或行就行不用反复重算。2.2 整体循环流程的伪代码主循环的思路大致是初始化结构模型参数设置地震动记录列表 设置强度水平向量 IM_target [0.05, 0.1, 0.2, ..., 1.5] for i 1:num_records 读取第i条地震动时程 标准化时程统一dt必要时截断 for j 1:length(IM_target) 计算该记录在当前目标强度下的缩放系数 缩放地震动 执行非线性时程分析 提取最大层间位移角 保存 IM_all(i,j) 和 EDP_all(i,j) end end这个两层循环是IDA计算的核心骨架。外层循环是地震动记录内层循环是强度水平。循环体内部的每一个环节都可以单独封装成函数比如scaleGroundMotion、runTimeHistoryAnalysis。这样做还有一个额外好处后期如果要并行化直接对第i条记录做parfor并行Matlab的并行计算工具箱就能派上用场。我在实际项目里一次IDA需要跑上百条工况如果每条工况的非线性时程分析用OpenSees之类的外部软件完成Matlab只负责生成输入文件、循环提交计算和读取结果那整个流程会慢一些但是批量处理能力很强。如果结构模型不太大直接在Matlab里做非线性时程分析省去了文件读写的时间速度反而更快。下面详细说每个关键环节的具体实现。3. 地震动缩放与强度指标选择看起来简单其实坑最多3.1 为什么优先用Sa(T1)而不是PGA很多新手习惯用PGA作为强度指标理由很充分PGA容易算、物理概念也清楚。但我在实际拟合易损性曲线时发现用PGA作为IM时不同记录之间的响应离散性非常大拟合出来的β对数标准差会偏大易损性曲线会显得特别平缓也就是“看起来什么都可能破坏什么都可能不破坏”。问题出在PGA只反应了地震动峰值的大小却没有反应地震动的频谱成分。一个长周期的框架结构如果输入一条高频成分很足但峰值很大的地震波结构可能根本没多大反应反过来一条峰值不大但含有较长周期成分的波反而可能激起结构强烈响应。应对方式是用结构基本周期T1处的谱加速度Sa(T1, 5%)作为IM。这相当于把地震动在“对结构最重要的频段”上的能量作为强度度量。这样做能显著降低记录间变异性让IDA曲线的“束”更集中拟合出的易损性曲线也更稳定。代价是你需要先算反应谱得到每条记录在T1处的谱值然后再进行缩放。3.2 两种缩放方法在Matlab里的实现按PGA缩放的计算公式很简单原始波最大绝对值如果是0.15g目标PGA是0.3g缩放系数就是0.3/0.152整条时程乘以2。注意这里的“g”是重力加速度在代码里通常把加速度单位统一成m/s²算完再换算。按Sa(T1)缩放稍微麻烦一点需要先计算原始记录在T1处的谱加速度。计算反应谱的常见办法是用一个小程序对单自由度体系做一系列不同周期下的时程响应取最大绝对加速度再乘上放大系数。Matlab里可以这么写一个简化版function [Sa, period] computeResponseSpectrum(acc, dt, zeta, Tlist) % acc : 加速度时程 (m/s^2) % dt : 时间步长 % zeta: 阻尼比 % Tlist: 需要计算谱值的周期列表 for i 1:length(Tlist) T Tlist(i); omega 2*pi/T; % 使用Newmark-beta法计算单自由度的位移响应 % 然后按位移响应与拟加速度的关系求Sa u newmarkSDOF(acc, dt, omega, zeta); Sa(i) max(abs(u)) * omega^2; end end得到原始记录的Sa0后缩放系数就是lambda targetSa / Sa0。实际操作中为了减少非线性分析时长很多人会先把地震动按单位PGA标准化再乘以目标PGA或目标Sa。代码层面就是两步pg_orig max(abs(acc)); acc_norm acc / pg_orig; % 单位PGA lambda targetPGA / 1.0; % 如果按PGA缩放 acc_scaled acc_norm * lambda; % 实际输入波一个必须提醒的点缩放同一个强度水平时每条记录缩放系数是完全不同的。有些记录的原始谱值在T1附近比较小缩放后峰值会大得离谱甚至超过几g这时候时程积分的收敛性和结构响应都可能失真。处理办法是适当调整强度水平的上限或者检查每条记录的缩放系数是否在一个合理范围内。3.3 时间步长和截断处理的细节地震动记录原始时间步长通常不统一有0.005s也有0.01s的。做非线性时程分析前必须先统一时间步长。Matlab里可以用resample函数但要注意重采样会引入滤波效应高频成分会受影响。我更建议的做法是在读取记录时就做线性插值统一成分析需要的dt比如0.005s或0.01s。另外很多记录后面有一段很长的“尾巴”加速度幅值已经很小了但对弹塑性结构来说尾部的地震波仍然可能让结构残余位移继续累积。要不要截断取决于你关心的响应是峰值还是残余值。如果只取最大层间位移角尾部作用不大但如果算残余位移就不能随便截断。我一般的做法是把加速度幅值小于峰值0.1%之后的段截掉但保留总时长的95%以上这样既控制计算时间也不影响峰值响应。4. 非线性时程分析模块的Matlab实现4.1 从单自由度弹塑性模型说起如果你只是想把IDA流程跑通第一步可以先做一个单自由度体系恢复力模型用理想弹塑性或者Bouc-Wen滞回模型。单自由度模型虽然结构简单但已经能完整走一遍“缩放→分析→提响应→累加→拟合”的全流程。等你把流程吃透再扩展成多自由度剪切模型或者接OpenSees外部模型思路是通用的。单自由度运动方程是这样的方程为m * ü c * u̇ F_s(u) -m * a_g(t)其中F_s是恢复力对理想弹塑性模型来说当位移小于屈服位移时是线性的k*u达到屈服力后保持恒定。这个方程用Newmark-β法可以稳定积分。下面是一段实现弹塑性单自由度时程分析的伪代码function [umax] nlSDOF(ag, dt, m, k, fy, zeta) % ag : 地震动加速度时程 % dt : 时间步长 % m : 质量 % k : 初始刚度 % fy : 屈服强度 % zeta: 阻尼比 omega sqrt(k/m); c 2 * m * omega * zeta; % Newmark参数平均加速度法 gamma 0.5; beta 0.25; u 0; v 0; fs 0; for i 1:length(ag)-1 % 预测-校正实现塑性状态的更新 % 先计算等效刚度与增量荷载 [u_new, v_new, fs_new] ... stepElasticPlastic(u, v, fs, ag(i), ag(i1), dt, m, c, k, fy, beta, gamma); u u_new; v v_new; fs fs_new; umax max(umax, abs(u)); end end这段代码只是骨架真正的stepElasticPlastic函数要判断当前位移是否超出屈服面超出后按塑性刚度更新恢复力需要用到状态变量并处理加卸载。对多自由度剪切模型每个楼层就是一个类似的单自由度单元但楼层之间有耦联方程变成矩阵形式。Matlab里组装质量、阻尼、刚度矩阵后用同样的Newmark方法直接解。4.2 多自由度剪切模型和Rayleigh阻尼的坑实际结构分析通常用多自由度模型起码3层以上才有意义。剪切模型是把每层质量集中到楼面处层间用弹塑性弹簧连接。刚度矩阵是经典的三对角形式。这个模型虽然简单但对评估框架结构的层间位移角已经足够用了。一个容易出问题的地方是阻尼。实际结构分析很少直接给质量比例阻尼而是用Rayleigh阻尼也就是阻尼矩阵写成质量和刚度的线性组合C αM βK。α和β的取值如果按前两阶阻尼比来标定就会导致一个非常不理想的结果——高阶振型的阻尼比可能被放大好几倍把高阶结构的真实响应压下去最大层间位移角会被低估。我在某次IDA计算里用固定瑞利阻尼参数跑了几个强度水平发现结构响应随强度增长几乎不增长检查半天才发现是阻尼在高阶振型上作用过强。后来改成只保留质量比例阻尼或者把β取小一点响应才正常。这个细节在Matlab里很难自动发现建议在跑大批量之前先取一条记录、两三个强度水平检查一下结构响应随IM增大的趋势是否符合预期再全量跑。4.3 如何提取稳定的EDP最大层间位移角时程分析完成后提取最大层间位移角不是简单地取每层位移差的最大绝对值。实际执行时要把位移时程序列都保留下来然后做层间位移差再对所有楼层、所有时间点取最大值。如果结构经历了较大的塑性变形还要考虑残余层间位移的影响但总层间位移角在IDA中一般就用弹塑性最大响应。% floor_disp(i, k) : 第i个时刻第k层的位移 EDP_all 0; for k 1:nFloor drift abs(floor_disp(:, k) - floor_disp(:, k-1)) / storyHeight; EDP_all max(EDP_all, max(drift)); end注意如果用的是层间位移角作为EDP不同极限状态阈值要按层间位移角来定义。比如中国规范里对框架结构常用1/550弹性和1/50弹塑性用这些阈值时要明确对应的破坏状态定义是否一致。不一致的阈值做对比没有意义。5. 从IDA曲线到易损性曲线数据清洗、极限状态定义与拟合5.1 极限状态的阈值怎么定极限状态可以理解成“破坏等级的门槛”。常见做法是定义四档LS1 轻微破坏、LS2 中等破坏、LS3 严重破坏、LS4 倒塌。对框架结构层间位移角阈值可以取1/200、1/100、1/50、1/25这类经验值但要结合结构类型酌情调整。阈值选择直接影响易损性曲线形状。阈值取得高曲线整体右移阈值取得低曲线左移。建议在代码里做成参数方便批量对比不同阈值下的结果。同时多个极限状态可以画在同一张图上形成易损性曲线族这比画单条曲线信息量大得多。5.2 点估法每个IM水平下的失效频率最简单的易损性统计方法是点估计法。对某个IM水平比如0.3g用该IM水平下所有地震动记录中的失效记录数除以总记录数就得到该点的失效概率估计。这个估计是非参数化的结果是一列散点。其实这是“在IM0.3g时做N次时程分析其中n次超过了极限状态阈值失效概率就是n/N”。如果N取20条记录那么失效概率估计的置信区间还是挺宽的尤其当n接近0或N时方差很大。所以点估法适合作初步判断最终还是需要拟合一条完整分布。5.3 对数正态分布的最大似然拟合工程上用的易损性曲线默认形式是假设极限状态对应的IM阈值服从对数正态分布。数学表达就是P(EDP ≥ LS | IM im) Φ( (ln(im) - ln(θ)) / β )这里的θ称为中位值也就是失效概率为50%时对应的IM值β是对数标准差反映离散程度。为什么要选对数正态分布因为IDA拟合出来的失效概率随IM增加通常呈现S形对数正态能很好地描述这种“低端趋近0、高端趋近1”的分布。最大似然估计的基本思路是在所有IM水平下已经知道每个IM点有N条记录、n条失效那么每点的观测结果服从二项分布把这些二项分布概率相乘就是似然函数。最大化这个函数的对数值就得到θ和β的估计。在Matlab里可以直接用fminsearchfunction [theta, beta] fitFragility(IM, n, N) % IM : 强度水平向量 % n : 每个IM下的失效记录数 % N : 每个IM下的总记录数 negLogLik (p) -( ... sum(n .* log(normcdf((log(IM) - p(1)) / p(2))) ... (N - n) .* log(1 - normcdf((log(IM) - p(1)) / p(2))) )); p0 [log(0.5), 0.4]; % 初始值 p_opt fminsearch(negLogLik, p0); theta exp(p_opt(1)); beta p_opt(2); end写这个函数时有几个坑需要说明。如果某个IM水平下全部记录都失效nN第二个对数项里就出现log(0)。解决方式有两种一是把n0或nN的点在似然函数中去掉因为它们的贡献在极限时有限二是对概率做一个微小的裁剪比如限制在[1e-6, 1-1e-6]。我一般用前一种简单直接。初始值的选择也很重要。θ初值可以选在失效概率约50%附近的IM值β初值取0.3~0.4通常都稳。如果初值偏差太大优化可能收敛到错误局部解拟合曲线会奇怪。担心这一点的话可以先用点估法扫一遍数据粗略看一下失效概率跨过50%的那个IM区间在哪然后把θ初值设在那里。5.4 绘制易损性曲线族的完整代码拟合出θ和β之后画曲线就是填空了。完整的一段绘图代码如下IM_plot linspace(0.01, 1.5, 200); colors lines(4); % 4个极限状态 labels {LS1, LS2, LS3, LS4}; figure; hold on; for i 1:4 theta theta_all(i, 1); beta beta_all(i, 2); P normcdf((log(IM_plot) - log(theta)) / beta); plot(IM_plot, P, LineWidth, 1.8, Color, colors(i, :)); end xlabel(IM 强度指标); ylabel(超越极限状态的概率); legend(labels, Location, east);如果数据量比较多还可以把点估法的散点也画出来拟合线一对比能直观看出最大似然拟合是否合理。手上项目里出现过拟合线在某些区间严重偏离散点的情况多半是β估计被极端点带跑偏了需要检查是不是某个IM点的失效记录数异常。6. 实操中容易踩的坑离散性大、拟合发散、记录数量不足6.1 地震动记录数量到底需要多少条易损性拟合对记录数量非常敏感。记录条数太少比如只有5条每个IM水平下的失效频率只有0、0.2、0.4这种稀疏取值最大似然估计也没有足够信息识别θ和β拟合结果会非常不稳定。我个人的经验是最少也要12条以上理想是20~30条。这个数量下每个IM水平才能算出比较平滑的失效频率曲线。但地震动记录数量增加的同时分析工况数量是线性增长的。假设20条记录、15个强度水平就是300次非线性时程分析。如果单次分析需要几十秒甚至几分钟总耗时会变得不可接受。解决办法一是并行计算二是减少强度水平数量但保留覆盖范围三是用更高效的简化模型。6.2 拟合参数发散时怎么调试使用fminsearch拟合的最大似然函数在数据分布极端的时候容易发散输出θ巨大无比或者β接近0。这种情况经常出现在低强度区结构响应几乎都低于阈值n/N在所有中低IM点都是0只有最高两个IM点开始有失效。此时数据携带的信息不足fminsearch就会往“β很小”的方向猛冲。我推荐的调试套路是两步走。第一步固定β为一个合理经验值比如0.35只拟合θ得到一个大致的θ。第二步把第一步得到的θ作为初值同时拟合θ和β。这样做能极大减少发散概率。另外如果拟合结果画出来后曲线在某个区间“跳崖”先怀疑是不是数据中某个IM点有异常记录回去检查该工况的时程分析是否收敛了。6.3 IM选择如何影响曲线形态用PGA和用Sa(T1)拟合出来的易损性曲线会有显著差异。通常用Sa(T1)时β更小曲线更陡中位值更集中用PGA时β更大曲线更平缓。这不是代码问题而是物理本质的体现。所以报告或者论文里一定要写清楚用的是哪个IM横轴单位是什么否则别人没法复现和对比。如果一定要用PGA做IM建议对结果多做一步敏感性分析看看β对记录集合变异的敏感程度。我在某次对比中发现仅更换5条记录用PGA时的中位值能差出15%以上而用Sa(T1)时只差3%到5%。这直接让我下定决心在后续项目里统一用Sa(T1)。6.4 对“全失效”和“全安全”强度水平的处理实际数据里最低的IM下所有记录都安全最高的IM下所有记录都失效这是正常现象。但这两个区域的点对拟合贡献和中间点不一样。中间区域的失效概率在0.1到0.9之间对θ和β识别最有效。全部安全的低强度点其实告诉算法“失效点在右侧”全部失效的高强度点告诉算法“失效点在左侧”它们的作用更像约束极值而不是精确插值。因此处理时不建议把所有IM点都塞进似然函数里而是保留那些至少有一个失效、又不全是失效的IM点。同时可以把每个IM点对应的记录数量N也作为加权信息写入似然函数这在方面上利用所有工况而不浪费。代码中可以在汇总阶段用一个逻辑索引过滤n0 nN的点再执行拟合过程。6.5 时间开销和批处理经验最后说一个和代码本身无关但决定项目体验的问题大批量IDA的时间开销。一次IDA跑300个工况如果每个工况用OpenSees耗时10秒那就是50分钟。但实际非线性分析经常不收敛调试起来会成倍放大时间成本。所以我的建议是先跑2条记录、3个强度水平快速验证代码流程再跑全部记录下的一两个中等强度水平检查分析收敛性最后全量跑。Matlab里可以用parfor替代for来并行外层循环。前提是每个循环迭代之间不共享复杂变量只写出结果。这里要注意Matlab的并行池在传递大数据时会有额外开销如果单条记录分析量很大建议在循环体内部尽量少保留中间变量只保留EDP结果内存和通讯压力都会小很多。我在实际项目里遇到过一次很有意思的情况某条记录在原强度水平下结构还处于弹性但稍微放大1.2倍后最大层间位移角突然跳了两倍多曲线出现一个“台阶”。这不是代码bug而是该记录恰好激活了结构的某个高阶响应分量或者进入了明显的塑性重分布阶段。这类现象在IDA曲线族里很常见不要急着删点先确认该工况的位移时程是否合理再决定是否保留。易损性曲线拟合本身对这种局部跳点不敏感但如果跳点出现在失效概率50%附近就需要额外认真检查。编写IDA代码的过程本质上就是不断在“计算效率和结果的合理性”之间找平衡。先把单自由度流程跑通再把模型复杂化每一步都留好检查输出最后得到的那张易损性曲线族图一定比直接套现成函数靠谱得多。
