尧图网页设计 
简介这套课程设计基于Python实现基本的地震易损性分析集成了完整源代码与配套文档面向土木工程、计算机、人工智能等专业的学生、教师及工程技术人员既能用于课设与毕设也适合作为Python工程分析的入门实战项目。项目压缩包共147个文件、约11.26MB其中核心代码为6个Python脚本16个eps与16个png图片呈现回归、残差及易损性曲线等分析结果9个txt文档提供说明与备注整体层次分明便于查找。目前已有80人学习浏览。代码均在测试运行成功后上传曾作为毕业设计答辩项目并获答辩评审平均分96分除源码外还提供了清晰的文档说明可帮助理解地震易损性分析从数据处理、回归拟合到曲线绘制的完整流程。在此基础上可自行修改扩展实现新的分析功能也能作为课程设计或项目立项演示的可靠基础。1. 为什么地震易损性分析要先看 ln 回归残差而不是直接画曲线拿到一个地震易损性分析的 Python 课程设计打开源码目录大概率会看到这样的输出文件ln(Pier1_cdrift) Regression.eps、ln(Bearing6_bdisp) Residual.eps、Abutment2_adispp FragilityCurve.eps。从文件名就能看出三层关系先把构件位移或延性取对数做最小二乘回归再检查残差是否随机最后把回归参数带进对数正态累计分布函数才能画易损性曲线。这里的反直觉点在于易损性曲线并不是“统计某次地震里多少结构损坏了”画出来的而是基于一个简化的需求模型外推得到ln(需求) a b * ln(地震动强度)。源代码里大量使用.eps矢量图说明作者的目标不是只跑通程序而是要把回归证据写进课程报告。下面直接盯住 ln 回归、残差检验、易损性曲线生成三个动作适合刚完成 Python 基础语法、准备把结构易损性分析落地的学生也适合想快速复现完整流程的工程师。2. 双参数对数正态模型易损性函数背后的回归依据2.1 为什么需求预测要选幂函数而不是线性函数桥梁构件的地震需求比如墩柱曲率延性Pier1_cdrift、桥台位移Abutment2_adispp、支座位移Bearing6_bdisp在低强度地震动下增长缓慢在强震下进入非线性后响应迅速放大。直接做线性回归时大震样本的离散度压过了小震样本回归线会被大震点“带偏”。用一个幂函数表达式描述需求与地震动强度指标IM的关系D a * IM^b * ε两边取对数后变成ln(D) ln(a) b * ln(IM) ln(ε)如果ln(ε)服从正态分布那么在ln(IM)-ln(D)空间里做最小二乘回归残差就是对数空间里的预测误差。这样处理有两个直接好处一是小震和大震的误差权重被拉平残差更容易满足等方差假设二是回归斜率 b 直接反映了构件需求对地震动强度的敏感性桥墩、桥台、支座之间可以横向比较。源代码中文件名前面的ln(...)前缀正是为了让读者意识到这张图的对数坐标系而不是文件名随便起的。2.2 从回归残差中提取标准差 β用scipy.stats.linregress拟合直线能拿到斜率、截距、R² 和 p 值但它不会直接返回残差标准差。残差标准差是后面易损性曲线斜率的关键参数必须自己计算。设拟合值为ln(D_i_pred) ln(a) b * ln(IM_i)残差为r_i ln(D_i) - ln(D_i_pred)所有残差的样本标准差记为 β_d代表需求对数空间的随机不确定性。在结构工程文献里这个 β_d 经常被直接叫做“回归标准差”或“对数标准差”。它决定了易损性曲线是陡峭还是平缓β_d 越小需求预测越精准脆弱性曲线在特定能力阈值下就越“干脆”β_d 越大表示同样地震强度下响应可能差别很大曲线斜率变缓。除了需求侧的不确定性构件能力本身也有离散性。工程上常用一个附加对数标准差 β_c 来描述比如钢筋混凝土墩柱的曲率延性能力离散、桥台位移限值的模型误差。最终计算失效概率时总标准差为β_tot sqrt(β_d^2 β_c^2)在课程设计中若没有能力试验数据β_c 常被设为 0.2~0.4。源代码里如果直接只给了 Regression 图和 Residual 图没有给能力标准差一般按 0.3 处理是比较常见的选择。2.3 能力阈值与失效概率公式的对数正态表达构件失效条件定义为“需求大于能力”。当需求和能力都是对数正态随机变量时给定 IM 下的失效概率可以写成标准正态累计分布函数P_fail(IM) Φ( (ln(a) b * ln(IM) - ln(C_m)) / β_tot )其中C_m是构件能力的中位值。实际项目里不同构件的C_m取值差距很大。下表列出三个典型构件的取值示例具体数值应以源代码里的说明为准构件输出名响应量能力中位示例说明Pier1_cdrift墩柱曲率延性比0.04延性限值查桥墩塑性铰模型Abutment2_adispp桥台位移m0.05桥台挡块或桩基容许位移Bearing6_bdisp支座位移m0.10常见板式橡胶支座剪切变形限值注意如果C_m的单位和D不一致比如位移取毫米、延性取无量纲比值回归出来的ln(a)会整体偏移易损性曲线中位值也会出错。起手的第一步应该是检查所有响应列的量纲统一而不是先跑回归。3. Python 实现从时程响应数据到回归图、残差图和易损性曲线3.1 数据结构与 CSV 组织源代码没有给出原始时程数据但一般课程设计中的数据表是csv或xlsx格式按“每行一次地震动输入每列一个响应量”排列。常见做法是这样工况IMPier1_cdriftAbutment2_adisppBearing6_bdispGM0010.120.00320.00510.0084GM0020.250.00810.01040.0162GM0030.510.02470.02330.0490第一列是地震波编号第二列是地震动强度指标IM后面各列是目标构件的最大响应。如果需要做多个指标的对比分析可以额外增加一列Sa_T1或PGA。读取时用 pandas 沿线解析即可。需要注意某些工况可能在中震级别下构件没屈服响应接近零取对数会出现-inf或者很夸张的负值。应对策略是对小于某极小阈值的数据保留但标记一般不要直接删掉因为这会改变回归样本的分布。3.2 用 scipy 做对数回归并计算残差标准差下面这段代码完成核心的回归计算。为了便于理解拆成数据读取、单构件拟合、结果输出三步import numpy as np import pandas as pd from scipy import stats # 读取时程响应记录 data pd.read_csv(seismic_im_d.csv) im data[IM].to_numpy() # 选择单个构件做回归分析 y data[Pier1_cdrift].to_numpy() # 取对数ln(IM) 与 ln(需求) x np.log(im) ylog np.log(y) # 最小二乘线性拟合 reg stats.linregress(x, ylog) # 计算残差和回归标准差 y_pred reg.intercept reg.slope * x residual ylog - y_pred beta np.std(residual, ddof2) print(fslope(b) {reg.slope:.4f}) print(fintercept(ln a) {reg.intercept:.4f}) print(fR^2 {reg.rvalue**2:.4f}) print(fresidual std beta_d {beta:.4f})逻辑说明linregress接受两个等长数组返回斜率slope、截距intercept、相关系数rvalue等。把IM和Pier1_cdrift分别取np.log等价于在双对数空间做线性拟合。ddof2表示计算残差标准差时减去两个自由度这是小样本下对回归模型自由度损失的一种修正课程报告里写这个细节会让答辩印象分提升。参数说明x是ln(IM)y是ln(Pier1_cdrift)。斜率b越接近 1.5说明墩柱延性需求随着地震动强度增长而急剧放大截距ln(a)是回归线在对数空间中的基准位置受IM单位影响很大。这里的启发在于若把IM从g改成cm/s²截距会变化但斜率不变因此跨项目对比时一定要写清楚强度指标单位。3.3 批量绘制 Regression 与 Residual 图并输出 eps课程设计里通常要处理多个构件逐个手写重复代码太低效。下面用一个列表循环生成所有构件的回归图和残差图。代码里使用 Matplotlib 的savefig(..., formateps)直接输出与源代码同名的.eps文件。import matplotlib.pyplot as plt components [Pier1_cdrift, Abutment2_adispp, Bearing6_bdisp] for comp in components: ylog np.log(data[comp].to_numpy()) reg stats.linregress(x, ylog) y_pred reg.intercept reg.slope * x residual ylog - y_pred # 两个子图左回归直线右残差分布 fig, (ax_left, ax_right) plt.subplots(1, 2, figsize(9, 4)) ax_left.scatter(x, ylog, s14, alpha0.6, labelobserved) x_dense np.linspace(x.min(), x.max(), 200) ax_left.plot(x_dense, reg.intercept reg.slope * x_dense, --, colorred, labelln regression) ax_left.set_xlabel(ln(IM)) ax_left.set_ylabel(fln({comp})) ax_left.legend() ax_right.scatter(x, residual, s14, alpha0.6) ax_right.axhline(0.0, colorblack, linewidth0.8) ax_right.set_xlabel(ln(IM)) ax_right.set_ylabel(residual) plt.tight_layout() plt.savefig(fln({comp})_Regression.eps, formateps) plt.close(fig)这段代码把每一个构件的回归结果和残差图写到同一个figure里导出为ln(Pier1_cdrift)_Regression.eps这种命名风格。x_dense用于画平滑的回归直线避免只画一条连接两个端点的生硬线段。残差图里如果看到点分布不是水平的带状而是开口喇叭形说明对数变换还不够可能需要考虑给IM换指数或分段拟合。3.4 易损性曲线绘图回归参数有了之后画脆弱性曲线的核心就是把P_fail(IM)公式转成代码。以下代码给出对每个构件分别绘制的实现from scipy.stats import norm cap_map { Pier1_cdrift: 0.04, Abutment2_adispp: 0.05, Bearing6_bdisp: 0.10, } beta_c 0.30 IM_range np.geomspace(im.min(), im.max(), 300) ln_IM np.log(IM_range) plt.figure(figsize(7, 5)) for comp in components: ylog np.log(data[comp].to_numpy()) reg stats.linregress(x, ylog) y_pred reg.intercept reg.slope * x residual ylog - y_pred beta_d np.std(residual, ddof2) beta_total np.sqrt(beta_d**2 beta_c**2) mu_lnD reg.intercept reg.slope * ln_IM capability np.log(cap_map[comp]) pf norm.cdf((mu_lnD - capability) / beta_total) plt.plot(IM_range, pf, lw2, labelcomp) plt.xscale(log) plt.xlabel(IM (g)) plt.ylabel(P(fail | IM)) plt.legend() plt.grid(True, linestyle--, linewidth0.4) plt.tight_layout() plt.savefig(FragilityCurve_All.eps, formateps) plt.savefig(FragilityCurve_All.png, dpi150)参数说明np.geomspace在对数坐标上生成 300 个 IM 点避免线性空间在小震区域取点太少。beta_total是需求残差标准差和能力标准差的平方和开方。norm.cdf的计算符号需要注意漏掉负号会导致曲线递增变成递减正确写法是需求超过能力时失效即(mu_lnD - ln(C_m)) / beta_total增大则失效概率增大。如果最终曲线随 IM 增大反而下降优先检查这个公式的正负号。cap_map里的能力中位值不对时整条曲线会在 IM 轴上左右平移但不改变曲线形状。课程设计中如果直接复制网上的阈值一定要结合源代码里的构件定义。比如某个桥墩已经算出塑性铰极限转角就不能再沿用默认延性比 0.04。4. 回归参数与 FragilityCurve 输出文件的对应关系、常见坑点排查4.1 文件名中的 Regression、Residual、FragilityCurve 对应什么源代码输出的.eps文件命名非常规则看文件名就能判断分析状态。下面用表格整理实际映射关系输出文件横轴纵轴告诉你什么ln(Abutment2_adispp) Regression.epsln(IM)ln(桥台位移)对数线性模型拟合是否成立ln(Abutment2_adispp) Residual.epsln(IM)回归残差残差是否存在趋势或异方差Abutment2_adispp FragilityCurve.epsIM失效概率构件易损性的最终判定结果ln(Pier1_cdrift) Regression.epsln(IM)ln(曲率延性)墩柱响应放大速度Pier1_cdrift FragilityCurve.epsIM失效概率墩柱倒塌或损伤曲线如果某个构件只出现Regression.eps而没有FragilityCurve.eps通常是因为缺少能力阈值C_m或者该构件的响应数据存在大量零值无法取对数。此时回到数据清洗把零值替换成极小正数如1e-5再试比直接把整条工况删掉要合理。4.2 回归系数 a、b 和 β 的解读方式以桥墩Pier1_cdrift为例如果拟合出的结果是这样一组示意值构件bln(a)β_dR²Pier1_cdrift1.48-2.960.320.86Abutment2_adispp1.06-3.710.250.74Bearing6_bdisp0.92-4.590.480.63b1.48说明墩柱延性需求对地震动强度很敏感ln(a)-2.96对应的是基准响应水平Bearing6_bdisp的 β_d 达到 0.48说明支座位移在相同地震强度下离散度更大可能是摩擦滑移或挡块碰撞导致。这种情况下再去要求易损性曲线非常陡峭是不现实的。桥梁抗震易损性分析里常见的误读是把 R² 当唯一质量标准。实际上R²0.74的桥台位移模型也可以继续使用因为需求预测模型不一定要求完全线性只要残差分布无偏、β_d 稳定即可。反过来如果 R² 很高但残差明显弯曲回归线只是强制穿过数据中间易损性曲线中段反而可能偏差很大。可以增加ln(IM)的平方项做多项式对比观察残差是否有明显系统性弯曲这是判断是否需要更高阶模型的简单办法。4.3 阈值与 β_c 调整影响同一个构件的能力中位值和能力标准差不同易损性曲线差异很大。以Bearing6_bdisp为例当能力中位C_m从 0.10 m 降到 0.06 m 时中位失效概率对应的 IM 会明显左移当 β_c 从 0.3 提高到 0.6曲线斜率明显放缓。课程报告的结论部分应该写明“为什么采用该阈值”常见依据是《公路桥梁抗震设计规范》中的支座容许位移或者文献里同类桥型的能力模型。源代码如果没有给出阈值来源建议在 README 或答辩幻灯片里补一句“能力参考自某规范/论文”这比列出公式更容易让老师认可。这里的实操检查方式是固定一组阈值只调整 β_c画出两组易损性曲线并观察差异。如果 β_c 从 0.2 变到 0.5曲线变化仍然很小说明需求侧不确定性占主导后续优化重点应放在提高回归质量上如果曲线变化剧烈说明能力不确定度不可忽略需要从材料本构或试验数据中取得更可靠的C_m。4.4 运行环境与常见报错在下载源代码后最容易卡住的三个问题eps文件打开为空白、中文乱码、scipy版本不兼容。.eps文件在 Windows 自带的图片查看器可能打不开推荐用 Acrobat Reader 或 Inkscape 查看如果是投稿论文直接嵌入 LaTeX 更合适。matplotlib 保存.eps时默认字体是 Type 3部分期刊会要求 Type 42 或 TrueType可以用plt.rcParams[pdf.fonttype] 42和plt.rcParams[ps.fonttype] 42规避。另一个常见报错是numpy新版本对np.float的移除导致旧代码崩溃需要把源码里的np.float替换为float。Python 环境搭建建议直接使用 Anaconda安装代码为pip install pandas scipy matplotlib或conda install ...可以顺手省去大量编译时间。5. 批量构件易损性分析与小样本置信区间计算技巧易损性曲线代码真正能继续深入的地方是把“单条曲线”升级为“带置信区间的曲线”。由于地震波样本数通常只有 20~50 条回归斜率、截距和残差标准差都有抽样误差。默认曲线是一条点估计不能回答“到底有多可信”这个问题。比较轻量的是用 Bootstrap 重采样来评估中位值 IM_50 的不确定范围。IM_50 是失效概率 0.5 对应的地震动强度用回归参数可以显式写出IM_50 exp( (ln(C_m) - ln(a)) / b )这个公式很简单却经常被忽略。它把易损性曲线的中位位置和回归参数直接挂钩只要有回归系数和能力阈值IM_50 就可以先算出来作为验证参照。如果从图里读出的 50% 概率点与公式算出的对不上一定是在norm.cdf里正负号或者beta_total合并出了问题。再看 Bootstrap 实现。对整个数据行做有放回抽样每次抽样拟合一次回归参数得到一条 IM_50 或指定 IM 下的失效概率重复 500~1000 次就能画出 5%~95% 置信带。核心代码如下rng np.random.default_rng(42) n_bootstrap 500 pf_matrix np.zeros((n_bootstrap, len(IM_range))) for i in range(n_bootstrap): idx rng.integers(0, len(data), sizelen(data)) sample data.iloc[idx] im_sample np.log(sample[IM].to_numpy()) y_sample np.log(sample[Bearing6_bdisp].to_numpy()) reg_i stats.linregress(im_sample, y_sample) residual_i y_sample - (reg_i.intercept reg_i.slope * im_sample) beta_d_i np.std(residual_i, ddof2) beta_total_i np.sqrt(beta_d_i**2 beta_c**2) mu_i reg_i.intercept reg_i.slope * ln_IM pf_matrix[i] norm.cdf((mu_i - np.log(cap_map[Bearing6_bdisp])) / beta_total_i) pf_lower np.percentile(pf_matrix, 5, axis0) pf_upper np.percentile(pf_matrix, 95, axis0)逻辑说明rng.integers生成随机行号有放回抽样后重复拟合。pf_matrix每一行是一次重采样得到的完整易损性曲线最后用percentile取出包络带。为什么不用直接对公式误差传递解析推导因为b和β_d之间并不独立Bootstrap 能够自然处理这种统计相关性。500 次重采样已经够课程设计展示趋势1000 次以上更平滑。这里有一个实际经验如果重采样后的置信区间在 50% 失效概率处宽度超过一倍 IM 区间说明地震波数量严重不足需要扩展谱匹配的地震波库而不是继续调节 β_c。反过来如果置信带很窄但回归残差图明显歪斜说明样本数量够但模型形式有问题此时应对比ln(IM)与IM两种自变量下的残差标准差选择更稳定的一种。这个做法可以同时应付课程报告里的“不确定度分析”问题和答辩老师的提问比单纯贴一张主曲线更有说服力。本文还有配套的精品资源点击获取 Python地震易损性分析：从对数回归残差到易损性曲线
简介这套课程设计基于Python实现基本的地震易损性分析集成了完整源代码与配套文档面向土木工程、计算机、人工智能等专业的学生、教师及工程技术人员既能用于课设与毕设也适合作为Python工程分析的入门实战项目。项目压缩包共147个文件、约11.26MB其中核心代码为6个Python脚本16个eps与16个png图片呈现回归、残差及易损性曲线等分析结果9个txt文档提供说明与备注整体层次分明便于查找。目前已有80人学习浏览。代码均在测试运行成功后上传曾作为毕业设计答辩项目并获答辩评审平均分96分除源码外还提供了清晰的文档说明可帮助理解地震易损性分析从数据处理、回归拟合到曲线绘制的完整流程。在此基础上可自行修改扩展实现新的分析功能也能作为课程设计或项目立项演示的可靠基础。1. 为什么地震易损性分析要先看 ln 回归残差而不是直接画曲线拿到一个地震易损性分析的 Python 课程设计打开源码目录大概率会看到这样的输出文件ln(Pier1_cdrift) Regression.eps、ln(Bearing6_bdisp) Residual.eps、Abutment2_adispp FragilityCurve.eps。从文件名就能看出三层关系先把构件位移或延性取对数做最小二乘回归再检查残差是否随机最后把回归参数带进对数正态累计分布函数才能画易损性曲线。这里的反直觉点在于易损性曲线并不是“统计某次地震里多少结构损坏了”画出来的而是基于一个简化的需求模型外推得到ln(需求) a b * ln(地震动强度)。源代码里大量使用.eps矢量图说明作者的目标不是只跑通程序而是要把回归证据写进课程报告。下面直接盯住 ln 回归、残差检验、易损性曲线生成三个动作适合刚完成 Python 基础语法、准备把结构易损性分析落地的学生也适合想快速复现完整流程的工程师。2. 双参数对数正态模型易损性函数背后的回归依据2.1 为什么需求预测要选幂函数而不是线性函数桥梁构件的地震需求比如墩柱曲率延性Pier1_cdrift、桥台位移Abutment2_adispp、支座位移Bearing6_bdisp在低强度地震动下增长缓慢在强震下进入非线性后响应迅速放大。直接做线性回归时大震样本的离散度压过了小震样本回归线会被大震点“带偏”。用一个幂函数表达式描述需求与地震动强度指标IM的关系D a * IM^b * ε两边取对数后变成ln(D) ln(a) b * ln(IM) ln(ε)如果ln(ε)服从正态分布那么在ln(IM)-ln(D)空间里做最小二乘回归残差就是对数空间里的预测误差。这样处理有两个直接好处一是小震和大震的误差权重被拉平残差更容易满足等方差假设二是回归斜率 b 直接反映了构件需求对地震动强度的敏感性桥墩、桥台、支座之间可以横向比较。源代码中文件名前面的ln(...)前缀正是为了让读者意识到这张图的对数坐标系而不是文件名随便起的。2.2 从回归残差中提取标准差 β用scipy.stats.linregress拟合直线能拿到斜率、截距、R² 和 p 值但它不会直接返回残差标准差。残差标准差是后面易损性曲线斜率的关键参数必须自己计算。设拟合值为ln(D_i_pred) ln(a) b * ln(IM_i)残差为r_i ln(D_i) - ln(D_i_pred)所有残差的样本标准差记为 β_d代表需求对数空间的随机不确定性。在结构工程文献里这个 β_d 经常被直接叫做“回归标准差”或“对数标准差”。它决定了易损性曲线是陡峭还是平缓β_d 越小需求预测越精准脆弱性曲线在特定能力阈值下就越“干脆”β_d 越大表示同样地震强度下响应可能差别很大曲线斜率变缓。除了需求侧的不确定性构件能力本身也有离散性。工程上常用一个附加对数标准差 β_c 来描述比如钢筋混凝土墩柱的曲率延性能力离散、桥台位移限值的模型误差。最终计算失效概率时总标准差为β_tot sqrt(β_d^2 β_c^2)在课程设计中若没有能力试验数据β_c 常被设为 0.2~0.4。源代码里如果直接只给了 Regression 图和 Residual 图没有给能力标准差一般按 0.3 处理是比较常见的选择。2.3 能力阈值与失效概率公式的对数正态表达构件失效条件定义为“需求大于能力”。当需求和能力都是对数正态随机变量时给定 IM 下的失效概率可以写成标准正态累计分布函数P_fail(IM) Φ( (ln(a) b * ln(IM) - ln(C_m)) / β_tot )其中C_m是构件能力的中位值。实际项目里不同构件的C_m取值差距很大。下表列出三个典型构件的取值示例具体数值应以源代码里的说明为准构件输出名响应量能力中位示例说明Pier1_cdrift墩柱曲率延性比0.04延性限值查桥墩塑性铰模型Abutment2_adispp桥台位移m0.05桥台挡块或桩基容许位移Bearing6_bdisp支座位移m0.10常见板式橡胶支座剪切变形限值注意如果C_m的单位和D不一致比如位移取毫米、延性取无量纲比值回归出来的ln(a)会整体偏移易损性曲线中位值也会出错。起手的第一步应该是检查所有响应列的量纲统一而不是先跑回归。3. Python 实现从时程响应数据到回归图、残差图和易损性曲线3.1 数据结构与 CSV 组织源代码没有给出原始时程数据但一般课程设计中的数据表是csv或xlsx格式按“每行一次地震动输入每列一个响应量”排列。常见做法是这样工况IMPier1_cdriftAbutment2_adisppBearing6_bdispGM0010.120.00320.00510.0084GM0020.250.00810.01040.0162GM0030.510.02470.02330.0490第一列是地震波编号第二列是地震动强度指标IM后面各列是目标构件的最大响应。如果需要做多个指标的对比分析可以额外增加一列Sa_T1或PGA。读取时用 pandas 沿线解析即可。需要注意某些工况可能在中震级别下构件没屈服响应接近零取对数会出现-inf或者很夸张的负值。应对策略是对小于某极小阈值的数据保留但标记一般不要直接删掉因为这会改变回归样本的分布。3.2 用 scipy 做对数回归并计算残差标准差下面这段代码完成核心的回归计算。为了便于理解拆成数据读取、单构件拟合、结果输出三步import numpy as np import pandas as pd from scipy import stats # 读取时程响应记录 data pd.read_csv(seismic_im_d.csv) im data[IM].to_numpy() # 选择单个构件做回归分析 y data[Pier1_cdrift].to_numpy() # 取对数ln(IM) 与 ln(需求) x np.log(im) ylog np.log(y) # 最小二乘线性拟合 reg stats.linregress(x, ylog) # 计算残差和回归标准差 y_pred reg.intercept reg.slope * x residual ylog - y_pred beta np.std(residual, ddof2) print(fslope(b) {reg.slope:.4f}) print(fintercept(ln a) {reg.intercept:.4f}) print(fR^2 {reg.rvalue**2:.4f}) print(fresidual std beta_d {beta:.4f})逻辑说明linregress接受两个等长数组返回斜率slope、截距intercept、相关系数rvalue等。把IM和Pier1_cdrift分别取np.log等价于在双对数空间做线性拟合。ddof2表示计算残差标准差时减去两个自由度这是小样本下对回归模型自由度损失的一种修正课程报告里写这个细节会让答辩印象分提升。参数说明x是ln(IM)y是ln(Pier1_cdrift)。斜率b越接近 1.5说明墩柱延性需求随着地震动强度增长而急剧放大截距ln(a)是回归线在对数空间中的基准位置受IM单位影响很大。这里的启发在于若把IM从g改成cm/s²截距会变化但斜率不变因此跨项目对比时一定要写清楚强度指标单位。3.3 批量绘制 Regression 与 Residual 图并输出 eps课程设计里通常要处理多个构件逐个手写重复代码太低效。下面用一个列表循环生成所有构件的回归图和残差图。代码里使用 Matplotlib 的savefig(..., formateps)直接输出与源代码同名的.eps文件。import matplotlib.pyplot as plt components [Pier1_cdrift, Abutment2_adispp, Bearing6_bdisp] for comp in components: ylog np.log(data[comp].to_numpy()) reg stats.linregress(x, ylog) y_pred reg.intercept reg.slope * x residual ylog - y_pred # 两个子图左回归直线右残差分布 fig, (ax_left, ax_right) plt.subplots(1, 2, figsize(9, 4)) ax_left.scatter(x, ylog, s14, alpha0.6, labelobserved) x_dense np.linspace(x.min(), x.max(), 200) ax_left.plot(x_dense, reg.intercept reg.slope * x_dense, --, colorred, labelln regression) ax_left.set_xlabel(ln(IM)) ax_left.set_ylabel(fln({comp})) ax_left.legend() ax_right.scatter(x, residual, s14, alpha0.6) ax_right.axhline(0.0, colorblack, linewidth0.8) ax_right.set_xlabel(ln(IM)) ax_right.set_ylabel(residual) plt.tight_layout() plt.savefig(fln({comp})_Regression.eps, formateps) plt.close(fig)这段代码把每一个构件的回归结果和残差图写到同一个figure里导出为ln(Pier1_cdrift)_Regression.eps这种命名风格。x_dense用于画平滑的回归直线避免只画一条连接两个端点的生硬线段。残差图里如果看到点分布不是水平的带状而是开口喇叭形说明对数变换还不够可能需要考虑给IM换指数或分段拟合。3.4 易损性曲线绘图回归参数有了之后画脆弱性曲线的核心就是把P_fail(IM)公式转成代码。以下代码给出对每个构件分别绘制的实现from scipy.stats import norm cap_map { Pier1_cdrift: 0.04, Abutment2_adispp: 0.05, Bearing6_bdisp: 0.10, } beta_c 0.30 IM_range np.geomspace(im.min(), im.max(), 300) ln_IM np.log(IM_range) plt.figure(figsize(7, 5)) for comp in components: ylog np.log(data[comp].to_numpy()) reg stats.linregress(x, ylog) y_pred reg.intercept reg.slope * x residual ylog - y_pred beta_d np.std(residual, ddof2) beta_total np.sqrt(beta_d**2 beta_c**2) mu_lnD reg.intercept reg.slope * ln_IM capability np.log(cap_map[comp]) pf norm.cdf((mu_lnD - capability) / beta_total) plt.plot(IM_range, pf, lw2, labelcomp) plt.xscale(log) plt.xlabel(IM (g)) plt.ylabel(P(fail | IM)) plt.legend() plt.grid(True, linestyle--, linewidth0.4) plt.tight_layout() plt.savefig(FragilityCurve_All.eps, formateps) plt.savefig(FragilityCurve_All.png, dpi150)参数说明np.geomspace在对数坐标上生成 300 个 IM 点避免线性空间在小震区域取点太少。beta_total是需求残差标准差和能力标准差的平方和开方。norm.cdf的计算符号需要注意漏掉负号会导致曲线递增变成递减正确写法是需求超过能力时失效即(mu_lnD - ln(C_m)) / beta_total增大则失效概率增大。如果最终曲线随 IM 增大反而下降优先检查这个公式的正负号。cap_map里的能力中位值不对时整条曲线会在 IM 轴上左右平移但不改变曲线形状。课程设计中如果直接复制网上的阈值一定要结合源代码里的构件定义。比如某个桥墩已经算出塑性铰极限转角就不能再沿用默认延性比 0.04。4. 回归参数与 FragilityCurve 输出文件的对应关系、常见坑点排查4.1 文件名中的 Regression、Residual、FragilityCurve 对应什么源代码输出的.eps文件命名非常规则看文件名就能判断分析状态。下面用表格整理实际映射关系输出文件横轴纵轴告诉你什么ln(Abutment2_adispp) Regression.epsln(IM)ln(桥台位移)对数线性模型拟合是否成立ln(Abutment2_adispp) Residual.epsln(IM)回归残差残差是否存在趋势或异方差Abutment2_adispp FragilityCurve.epsIM失效概率构件易损性的最终判定结果ln(Pier1_cdrift) Regression.epsln(IM)ln(曲率延性)墩柱响应放大速度Pier1_cdrift FragilityCurve.epsIM失效概率墩柱倒塌或损伤曲线如果某个构件只出现Regression.eps而没有FragilityCurve.eps通常是因为缺少能力阈值C_m或者该构件的响应数据存在大量零值无法取对数。此时回到数据清洗把零值替换成极小正数如1e-5再试比直接把整条工况删掉要合理。4.2 回归系数 a、b 和 β 的解读方式以桥墩Pier1_cdrift为例如果拟合出的结果是这样一组示意值构件bln(a)β_dR²Pier1_cdrift1.48-2.960.320.86Abutment2_adispp1.06-3.710.250.74Bearing6_bdisp0.92-4.590.480.63b1.48说明墩柱延性需求对地震动强度很敏感ln(a)-2.96对应的是基准响应水平Bearing6_bdisp的 β_d 达到 0.48说明支座位移在相同地震强度下离散度更大可能是摩擦滑移或挡块碰撞导致。这种情况下再去要求易损性曲线非常陡峭是不现实的。桥梁抗震易损性分析里常见的误读是把 R² 当唯一质量标准。实际上R²0.74的桥台位移模型也可以继续使用因为需求预测模型不一定要求完全线性只要残差分布无偏、β_d 稳定即可。反过来如果 R² 很高但残差明显弯曲回归线只是强制穿过数据中间易损性曲线中段反而可能偏差很大。可以增加ln(IM)的平方项做多项式对比观察残差是否有明显系统性弯曲这是判断是否需要更高阶模型的简单办法。4.3 阈值与 β_c 调整影响同一个构件的能力中位值和能力标准差不同易损性曲线差异很大。以Bearing6_bdisp为例当能力中位C_m从 0.10 m 降到 0.06 m 时中位失效概率对应的 IM 会明显左移当 β_c 从 0.3 提高到 0.6曲线斜率明显放缓。课程报告的结论部分应该写明“为什么采用该阈值”常见依据是《公路桥梁抗震设计规范》中的支座容许位移或者文献里同类桥型的能力模型。源代码如果没有给出阈值来源建议在 README 或答辩幻灯片里补一句“能力参考自某规范/论文”这比列出公式更容易让老师认可。这里的实操检查方式是固定一组阈值只调整 β_c画出两组易损性曲线并观察差异。如果 β_c 从 0.2 变到 0.5曲线变化仍然很小说明需求侧不确定性占主导后续优化重点应放在提高回归质量上如果曲线变化剧烈说明能力不确定度不可忽略需要从材料本构或试验数据中取得更可靠的C_m。4.4 运行环境与常见报错在下载源代码后最容易卡住的三个问题eps文件打开为空白、中文乱码、scipy版本不兼容。.eps文件在 Windows 自带的图片查看器可能打不开推荐用 Acrobat Reader 或 Inkscape 查看如果是投稿论文直接嵌入 LaTeX 更合适。matplotlib 保存.eps时默认字体是 Type 3部分期刊会要求 Type 42 或 TrueType可以用plt.rcParams[pdf.fonttype] 42和plt.rcParams[ps.fonttype] 42规避。另一个常见报错是numpy新版本对np.float的移除导致旧代码崩溃需要把源码里的np.float替换为float。Python 环境搭建建议直接使用 Anaconda安装代码为pip install pandas scipy matplotlib或conda install ...可以顺手省去大量编译时间。5. 批量构件易损性分析与小样本置信区间计算技巧易损性曲线代码真正能继续深入的地方是把“单条曲线”升级为“带置信区间的曲线”。由于地震波样本数通常只有 20~50 条回归斜率、截距和残差标准差都有抽样误差。默认曲线是一条点估计不能回答“到底有多可信”这个问题。比较轻量的是用 Bootstrap 重采样来评估中位值 IM_50 的不确定范围。IM_50 是失效概率 0.5 对应的地震动强度用回归参数可以显式写出IM_50 exp( (ln(C_m) - ln(a)) / b )这个公式很简单却经常被忽略。它把易损性曲线的中位位置和回归参数直接挂钩只要有回归系数和能力阈值IM_50 就可以先算出来作为验证参照。如果从图里读出的 50% 概率点与公式算出的对不上一定是在norm.cdf里正负号或者beta_total合并出了问题。再看 Bootstrap 实现。对整个数据行做有放回抽样每次抽样拟合一次回归参数得到一条 IM_50 或指定 IM 下的失效概率重复 500~1000 次就能画出 5%~95% 置信带。核心代码如下rng np.random.default_rng(42) n_bootstrap 500 pf_matrix np.zeros((n_bootstrap, len(IM_range))) for i in range(n_bootstrap): idx rng.integers(0, len(data), sizelen(data)) sample data.iloc[idx] im_sample np.log(sample[IM].to_numpy()) y_sample np.log(sample[Bearing6_bdisp].to_numpy()) reg_i stats.linregress(im_sample, y_sample) residual_i y_sample - (reg_i.intercept reg_i.slope * im_sample) beta_d_i np.std(residual_i, ddof2) beta_total_i np.sqrt(beta_d_i**2 beta_c**2) mu_i reg_i.intercept reg_i.slope * ln_IM pf_matrix[i] norm.cdf((mu_i - np.log(cap_map[Bearing6_bdisp])) / beta_total_i) pf_lower np.percentile(pf_matrix, 5, axis0) pf_upper np.percentile(pf_matrix, 95, axis0)逻辑说明rng.integers生成随机行号有放回抽样后重复拟合。pf_matrix每一行是一次重采样得到的完整易损性曲线最后用percentile取出包络带。为什么不用直接对公式误差传递解析推导因为b和β_d之间并不独立Bootstrap 能够自然处理这种统计相关性。500 次重采样已经够课程设计展示趋势1000 次以上更平滑。这里有一个实际经验如果重采样后的置信区间在 50% 失效概率处宽度超过一倍 IM 区间说明地震波数量严重不足需要扩展谱匹配的地震波库而不是继续调节 β_c。反过来如果置信带很窄但回归残差图明显歪斜说明样本数量够但模型形式有问题此时应对比ln(IM)与IM两种自变量下的残差标准差选择更稳定的一种。这个做法可以同时应付课程报告里的“不确定度分析”问题和答辩老师的提问比单纯贴一张主曲线更有说服力。本文还有配套的精品资源点击获取