简介：这是一份基于Python实现的地震易损性分析源码包，面向土木工程、地震工程方向的研究生、科研人员及结构设计人员，解决从地震需求计算到易损性曲线绘制的代码实现难题。压缩包共147个文件，包含6个Python源程序、100个out结果数据文件、16个eps矢量图、16个png图片和9个txt文本说明，整体大小11.29MB。其中eps与png清晰呈现了ln(工程需求参数)回归、残差分布及不同构件的易损性曲线，txt文件则提供配套说明，便于对照运行。已有372人学习下载，代码附有详细注释，下载即可直接运行，适合快速上手。内容预览表明覆盖桥台位移、桥墩延性比、支座位移等关键构件的回归建模与易损性曲线生成，对桥梁地震易损性分析、地震风险评估类课题具有较高参考价值。
1. 地震易损性分析为什么值得用 Python 重写一遍带注释的源码
一栋 20 层钢筋混凝土框架，设防 8 度，当地表峰值加速度到 0.3g 时，它进入中等破坏的概率是多少？这个数不是拍脑袋来的——把几十条地震动输入结构做非线性时程分析，再把结果统计成易损性曲线，这就是地震易损性分析的主体。重型分析在 OpenSees、ETABS 里完成，数据清洗、损伤判定、曲线拟合、批量出图，则由 Python 源代码收尾。这样一份带注释的源码适合两类人：土木方向的研究生，想拿到能改参数直接跑的管线；做抗震风险评估或保险定价的工程师，需要把易损性参数整理成可供外部平台调用的结果。下面按模型、数据处理、拟合、验证的顺序，把源码里最容易写错也最值得看的部分拆开讲。
2. 易损性曲线的对数正态模型与 Python 数据结构
拿到这类压缩包，我建议先画数据流再动代码：地震动记录和 IDA 反应文件进来，中间做层间位移角提取与损伤状态判定，最后得到每条记录的容量 IM，再拟合成对数正态曲线。这一章先把数学假设和数据结构钉死，后面的代码才能逐段核对。
2.1 为什么默认假设易损性服从对数正态分布
结构反应与地震动强度在一阶近似下呈幂律关系：层间位移角 θ 与阻尼比无关地近似满足 θ ≈ a·(Sa(T1))^b，两边取对数后是线性回归形式。地震动记录之间的差异体现在回归残差里，残差近似正态，于是"结构达到某个损伤限值所需的地震动强度"（容量 IM）就近似服从对数正态分布。易损性曲线的标准写法是：
F(x) = Φ( (ln x − λ) / β )
其中 λ 是 ln(容量 IM) 的均值，β 是对数标准差，Φ 是标准正态 CDF。λ 和 β 决定了整条曲线：λ 控制中位数的位置，β 控制曲线斜率，β 越大曲线越平缓，说明不同地震动之间结构反应越离散。
拟合 λ、β 有两条路线。一条是逐记录提取超限 IM，再对这批数做分布拟合，参数少、评审容易，注释也好写；另一条是先对 ln θ 和 ln IM 做回归，再把回归残差和限值的不确定性合成进 β，适合写论文但工程平台用得少。我一般建议源码主线走第一条，第二条保留成可选模块，标题里这类"源码+注释"项目也通常按第一条组织。
2.2 IDA 结果文件的典型列结构
易损性分析不要求自己写有限元，常见做法是用 OpenSees 跑完增量动力分析（IDA），把结果导出成文本或 CSV，Python 只做下游处理。文件结构通常是三张表：
| 文件 | 关键字段 | 在源码里的用途 |
|---|---|---|
| ground_motion.csv | record_id, pga_g, sa_t1_g | 每条记录的谱加速度，决定 IDA 缩放档位 |
| ida_summary.csv | record_id, scale_factor, max_drift_rad | 单次分析的峰值层间位移角，主输入 |
| out_drift_*.txt | time, drift1, drift2, ... | recorder 原始输出，供复核和重算 |
有项目把后两张表合并成宽表，有些还把缩放档位当成新的记录条数，拟合出来的 β 虚高到 0.8 以上，这是源码注释里要专门提醒的一条；反过来，如果只保留每条记录最大档位的反应，容量 IM 的插值就没法做，曲线中位数也会偏大。
2.3 用 dataclass 把三个核心对象装起来
源码里我习惯先用 dataclass 定义数据形态，而不是到处传裸字典；字段名带单位、注释写明含义，拟合和画图模块就共用同一套对象，后面加字段（比如加隔震层反应）时也只用改一处：
from dataclasses import dataclass from pathlib import Path @dataclass class GroundMotion: """一条地震动记录：唯一标识与缩放后的强度指标""" record_id: str pga_g: float = 0.0 # 缩放档位对应的峰值加速度，单位 g sa_t1_g: float = 0.0 # 结构基本周期 T1 处谱加速度，单位 g source_file: Path | None = None # 原始时程文件路径，复核时用 @dataclass class IDARecord: """单条记录在某个缩放档位下的反应结果""" gm: GroundMotion scale_factor: float = 1.0 max_drift_rad: float = 0.0 # 该档位整个时程的最大层间位移角，弧度 converged: bool = True # 非线性分析是否收敛，False 要单独处理用 dataclass 的好处是可以用 dataclasses.asdict 直接转成 DataFrame 或 JSON，也方便在注释里写明每个字段的单位。max_drift_rad 把单位写进字段名，能避免后面把 rad 和百分数混用，这是易损性分析源码里最常见的隐性 bug 来源。
3. 地震易损性分析的 IDA 数据整理：层间位移角与损伤状态判定
易损性分析只关心每个缩放档位下结构反应的峰值，不关心峰值出现的时刻或地震动中段细节，所以数据处理的第一步是压缩：把整条时程压成一个最大层间位移角，再把它压成一个损伤状态。这一章的所有函数都只做提取与判定，不碰拟合，边界保持干净。
3.1 从 recorder 时程里提取最大层间位移角
OpenSees 的 recorder 输出通常是纯文本，第一列是时间，后面每列是一个楼层或一个单元的位移角。常见做法是先用 pandas 按空白分隔读取，跳过文件头，再对整个数值区取绝对值最大值：
import numpy as np import pandas as pd def extract_max_drift(filepath: str, skiprows: int = 8) -> float: """读取一条时程，返回最大层间位移角（rad）""" raw = pd.read_csv(filepath, sep=r"\s+", skiprows=skiprows, header=None) drift_data = raw.iloc[:, 1:] # 首列是时间 max_drift = float(np.max(np.abs(drift_data.to_numpy()))) return max_driftskiprows 要按 recorder 输出头的实际行数调，取少了会吞掉数据，取多了直接报错。np.abs 是必要的，层间位移角有正负，取最大绝对值才代表最不利时刻。如果一条文件只记录单层位移角，drift_data 只有一列，这段代码同样成立。用整条时程取峰值而不是用包络值，是避免把静力推覆分析的结果混进易损性样本。
提示：OpenSees 输出默认单位是弧度，损伤状态阈值表里的 0.002 也是弧度，两边一致才能直接比。常见错误是把输出当百分数，阈值按 0.5% 去写，整条曲线会系统性左移。
3.2 损伤状态阈值：一套可改的配置而不是硬编码
四状态划分（轻微、中等、严重、倒塌）在文献里有多种取值，差异主要来自结构类型和规范背景。下面是针对 RC 框架较常见的一组，可作为默认值：
| 损伤状态 | 层间位移角 (rad) | 约合表述 | 常见出处/用途 |
|---|---|---|---|
| slight | 0.002 | 1/500 | 开裂初判，接近弹性限值 |
| moderate | 0.005 | 1/200 | 可修复性界限 |
| extensive | 0.015 | 1/67 | 与 FEMA 356 的 LS 同量级 |
| collapse | 0.040 | 1/25 | 倒塌判定常用下限 |
中国规范体系下，通常把"小震不坏"的弹性层间位移角取 1/550（约 0.0018），把"大震不倒"的弹塑性限值取 1/50（0.02）。前者可以作为 slight 的参考下限，后者用于校核整体稀有水准，不适合直接替换 collapse 阈值。处理这类源码时，阈值要集中放配置文件并在注释里写明出处，因为阈值一动，拟合出的中位数和 β 全变：
DAMAGE_LIMITS = { "slight": 0.002, # 参考 RC 框架经验取值，出处写进注释 "moderate": 0.005, "extensive": 0.015, "collapse": 0.040, }3.3 判定损伤状态并提取每条记录的容量 IM
有了每条记录的峰值漂移，下一步做两件事：一是给它打上损伤状态标签，用于汇总各状态占比；二是按缩放档位找到首次达到限值时的 IM，作为该记录的容量观测值：
def assign_damage_state(drift: float) -> str: """按字典序逐级比较，返回不超过该限值的最高损伤档位""" for state, limit in DAMAGE_LIMITS.items(): if drift < limit: return state return "collapse" # 超过所有限值才落到倒塌 def capacity_im(sa_levels: np.ndarray, drift_levels: np.ndarray, limit: float) -> float: """在单条记录的 IDA 档位里，线性插值求首次达到限值的 IM""" idx = np.where(drift_levels >= limit)[0] if len(idx) == 0: return np.inf # 到最大档位仍未达到，右截尾 i = idx[0] if i == 0: return float(sa_levels[0]) # 从最小档位就超限，保守取该值 return float(np.interp(limit, drift_levels[i-1:i+1], sa_levels[i-1:i+1]))assign_damage_state 依赖字典插入序，所以 config 里不要用无序 dict 或在外层重新排序。capacity_im 要求档位严格升序，np.interp 也要求 xp 单调，如果 IDA 曲线在中段出现下降（结构软化后的回复现象），取首次穿越档位而不是末段。返回 np.inf 的记录不能直接丢弃，它是右截尾样本，会在第四章的似然函数里起作用，源码注释里必须把这个约定写清楚。
4. 用 Python 最大似然估计拟合易损性曲线参数并出图
前两章得到的是每条记录、每个损伤状态的容量 IM 集合，第四章把它们拟合成 λ、β，并画出带置信区间的易损性曲线。这一章是源码的数学核心，也是代码注释最密集的部分。
4.1 为什么最大似然估计比直接画经验点更合适
一条常见错误路线是：把记录按容量 IM 排序，画出超限频率和 IM 的散点，再用 polyfit 去配平滑线。这样得到的是经验点而非分布参数，无法外推，也无法给出参数不确定性。容量法下，容量 IM 本身是一组独立观测值，假设服从对数正态分布后，MLE 的闭式解非常干净：
λ̂ = mean(ln x_i)，β̂ = std(ln x_i, ddof=1)
闭式解的优点是样本量 10 到 100 时都稳定，注释里也好说明：对容量 IM 取对数，算均值与标准差即可。真正的麻烦来自截尾：如果 IDA 最大档位不够高，部分记录到顶都没达到限值，观测值是 np.inf，这时必须用截尾似然，把"未超限"的信息通过生存函数计入，否则中位数会被系统性低估。
4.2 无截尾与右截尾两套拟合的实现
import numpy as np from scipy import stats, optimize def mle_lognormal(caps: np.ndarray) -> dict: """无截尾时的闭式解拟合, 返回中位数与对数标准差""" caps = np.asarray(caps, dtype=float) valid = caps[np.isfinite(caps)] if len(valid) < 3: raise ValueError("有效容量样本少于 3 条，拒绝拟合") lam = np.log(valid).mean() beta = np.log(valid).std(ddof=1) return {"median": float(np.exp(lam)), "beta": float(beta)} def mle_lognormal_censored(caps: np.ndarray) -> dict: """含右截尾样本的 MLE：np.inf 表示到最大档位仍未达到限值""" caps = np.asarray(caps, dtype=float) observed = caps[np.isfinite(caps)] censored = caps[~np.isfinite(caps)] if len(censored) == 0: return mle_lognormal(caps) def neg_loglik(theta): lam, beta = theta if beta <= 0: return 1e12 ll = stats.norm.logpdf(np.log(observed), lam, beta).sum() ll += stats.norm.logsf(np.log(censored), lam, beta).sum() return -ll x0 = [float(np.log(observed).mean()), 0.4] res = optimize.minimize(neg_loglik, x0, method="Nelder-Mead") return {"median": float(np.exp(res.x[0])), "beta": float(res.x[1])}第一段代码里，np.isfinite 过滤掉 inf，闭式解直接对 ln(x) 求均值和无偏标准差，ddof=1 是为了小样本下 β 不被低估。第二段是带截尾的似然，censored 数组里的每个值代表"容量 IM 至少大于该档位"，因此用 logsf 加上对数生存概率，再把整个负对数似然交给 minimize。初始值 x0 用观测样本的 log 均值，β 起步取 0.4，Nelder-Mead 不依赖梯度，对这类平滑目标足够稳定。
提示：截尾占比超过 20% 时，曲线右尾基本靠外推，评审时要显式说明。遇到这种情况优先回去把 IDA 的最大缩放档位提高，而不是在拟合上反复调参。
4.3 Bootstrap 置信区间与易损性曲线绘图
单次拟合只给一组点估计，工程报告里需要知道中位数的区间。常见做法是对观测样本做 2000 次有放回重采样，每次都重新拟合，取 2.5% 和 97.5% 分位作为 95% 区间：
def bootstrap_ci(caps: np.ndarray, n_boot: int = 2000, seed: int = 42): """自助重采样得到 median 与 beta 的 95% 置信区间""" rng = np.random.default_rng(seed) caps = np.asarray(caps, dtype=float) caps = caps[np.isfinite(caps)] medians, betas = [], [] for _ in range(n_boot): sample = caps[rng.integers(0, len(caps), len(caps))] medians.append(float(np.exp(np.log(sample).mean()))) betas.append(float(np.log(sample).std(ddof=1))) return np.percentile(medians, [2.5, 97.5]), np.percentile(betas, [2.5, 97.5]) def plot_fragility(median: float, beta: float, im_grid: np.ndarray, ax): """画一条标准对数正态易损性曲线, im_grid 单位与 IM 一致""" frac = stats.lognorm(s=beta, scale=median).cdf(im_grid) ax.plot(im_grid, frac, lw=2) ax.set_xlabel("PGA / g") ax.set_ylabel("P(DS ≥ ds | IM)") return frac画完之后有三个自检：中位数处曲线应穿过 50% 上下；β 落在 0.2 到 0.6 之间，β 小于 0.2 多半是损伤阈值设太低导致所有记录都在低 IM 处集中超限；im_grid 的起点不要取 0，对数正态在 0 处无定义，工程上从 0.05g 起步更实际。
5. 验证这套 Python 地震易损性分析源码的拟合管线
源码发到评审手里之前，先跑三件小事：合成数据冒烟测试、三个数据坑的排查、以及把结果导出成平台能直接用的格式。整套源码依赖就 numpy、scipy、pandas、matplotlib 四个包，装好 Python 后 pip install 一行到位，不需要额外的环境配置。
5.1 用已知参数的合成数据做冒烟测试
from scipy import stats import numpy as np true_median, true_beta = 0.35, 0.45 caps = stats.lognorm(s=true_beta, scale=true_median).rvs(200, random_state=7) res = mle_lognormal(caps) assert abs(res["median"] - true_median) < 0.03 assert abs(res["beta"] - true_beta) < 0.05如果断言失败，优先查 mle_lognormal 的输入是否被 pandas 解析成整型、或者 dataclass 的字段名被 asdict 改名导致长度不匹配；这类问题在带注释的源码里最常见，示例数据本身往往没问题。
5.2 三个最容易把整条曲线带歪的坑
- 阈值与单位不匹配：recorder 输出的是弧度，阈值表写的也是弧度，中间只要有一步把其中之一当成了百分数，整条曲线就会整体平移，拟合出的中位数直接失真。
- 跨记录插值：容量 IM 必须在单条记录内部逐档位插值，把所有记录的所有档位混在一起插值，等于把记录间差异和档位差异混成同一个离散度，β 会虚高到 0.8 以上。
- 非收敛档位整条删除：IDA 最后一级不收敛往往说明结构已经接近局部失效，把记录整条删掉会让曲线右尾被截短，保守做法是把这个档位的漂移判为超限，让它在拟合里体现出来。
5.3 导出 JSON 给上层风险平台调用
源码的最终交付不只是一张 png 图，建议把拟合参数连同数据来源一起导出成 JSON，上层做灾害损失评估或保险定价的平台拿到后可以直接用：
import json payload = { "im_type": "PGA", "unit": "g", "damage_state": "moderate", "median_g": res["median"], "beta": res["beta"], "model": "lognormal", "n_records": int(len(caps)), "censored_ratio": float(round((caps == np.inf).mean(), 3)), } with open("fragility_moderate.json", "w", encoding="utf-8") as f: json.dump(payload, f, indent=2, ensure_ascii=False)censored_ratio 是评审最关心的元数据，它直接告诉下游这条曲线的右尾有多少外推成分；平台拿到 JSON 后直接调 scipy.stats.lognorm 的 cdf 就能生成任意强度下的破坏概率，不需要再跑一次 Python 拟合管线。
本文还有配套的精品资源，点击获取
