简介针对地震工程分析需求这份MATLAB程序基于《抗震工程学》理论编写面向科研人员与在校学生可快速由地震波原始数据生成相对速度谱、相对位移谱和绝对加速度谱。资源包内共2个文件m源码文件用于执行谱计算与绘图txt文件存放支持运行的地震波数据整体压缩包仅11KB轻量易用。目前已有463人学习下载程序调试完备、精度可靠能够满足地震分析、课程作业及入门实践等场景。完成下载后即可直接加载数据运行通过结果直观对比三类反应谱曲线帮助理解结构地震响应特性也为后续研究提供了可扩展的代码基础和真实数据支撑。1. 同一张反应谱图为什么要分相对位移、相对速度和绝对加速度三条线拿到一条采样间隔 0.005s 的强震加速度记录时结构工程师通常不先问 PGA 是多少而是问周期 1s 左右的框架这条波能把它推到多大力。直接看波形回答不了——同一条波对周期 0.2s 和 2s 的单自由度体系的放大效应可能差一个数量级。反应谱就是为这个场景设计的对一组周期逐次求解单自由度运动方程提取相对位移、相对速度、绝对加速度的时程峰值得到相对位移谱 Sd、相对速度谱 Sv、绝对加速度谱 Sa 三条谱曲线。这里给出的 Matlab 程序覆盖从运动方程、Newmark 积分、参数标定到画图验证的完整链路适合做抗震分析、结构设计复核或课程作业的读者直接参照落地。2. 从运动方程到三条谱曲线反应谱在算什么2.1 单自由度体系的相对运动方程要理解相对位移谱跟绝对加速度谱的差别关键是把“相对”和“绝对”两个坐标系分清。取结构底部跟着地面一起平移结构质量 m 相对地面的位移记作 u(t)地面自身位移记作 ug(t)那么质量块的绝对位移是 u(t)ug(t)。以相对位移为未知量单自由度体系的运动方程为$$m\ddot{u}(t)c\dot{u}(t)ku(t)-m\ddot{u}_g(t)$$右侧的 $\ddot{u}_g(t)$ 就是地震记录里的加速度时程可以直接从 csv 或文本文件读进来。它不是传统意义上的外荷载而是地面运动引起的等效惯性力。这里的阻尼 c 和刚度 k 都不是随意给定的刚度 k 对应结构周期$km(2\pi/T)^2$阻尼 c 对应阻尼比$c2m\xi\omega$。实际程序里通常把质量 m 归一化成 1只保留阻尼比 $\xi$ 和圆频率 $\omega$ 两个自由度。绝对加速度在这里不是 $\ddot{u}(t)$而是 $\ddot{u}(t)\ddot{u}_g(t)$也就是质点相对地面的加速度加上地面本身的加速度。很多初学者以为反应谱里的“绝对加速度谱”指的是相对加速度的峰值这会导致后期跟规范里的地震影响系数对不上需要格外留意。2.2 相对位移谱、相对速度谱与绝对加速度谱的定义程序要输出的三条谱在文献里有固定的英文缩写和计算方式Sd 是相对位移时程绝对值的最大值Sv 是相对速度时程绝对值的最大值Sa 是绝对加速度时程绝对值的最大值。对同一条地震波不同的周期 T 会得到不同的 Sd、Sv、Sa 值把所有这些值连成曲线就构成反应谱。工程上还有一个经常出现的计算变体就是伪速度谱和伪加速度谱分别定义为 $\omega S_d$ 和 $\omega^2 S_d$。教科书和规范里的速度反应谱多数是伪速度谱而不是真实相对速度谱。两者在小阻尼下数值非常接近但定义不同。本程序的默认输出是严格定义值$S_d\max_t|u(t)|$$S_v\max_t|\dot u(t)|$$S_a\max_t|\ddot u(t)\ddot u_g(t)|$同时把伪谱作为校验手段保留在验证环节。2.3 Newmark-β 平均加速度法的 Matlab 实现时域积分使用 Newmark-β 法里最常用的平均加速度格式$\gamma1/2$、$\beta1/4$。它无条件稳定可以使用原始采样步长而不必为了稳定性加密积分步长。下面是单自由度体系反应谱计算函数输入一条加速度时程、采样间隔、阻尼比和周期输出三个谱值。function [Sd, Sv, Sa] sdof_response_spectrum(ag, dt, xi, T) % ag : 地面加速度时程单位任意cm/s^2 或 m/s^2 % dt : 采样间隔单位 s % xi : 阻尼比例如 0.05 % T : 单自由度体系周期单位 s % Sd Sv Sa 分别为相对位移、相对速度、绝对加速度峰值 omega 2 * pi / T; % 圆频率 m 1; % 质量归一化 c 2 * xi * omega * m; % 阻尼系数 k omega^2 * m; % 刚度系数 % Newmark 平均加速度常数 gamma 1 / 2; beta 1 / 4; n length(ag); u zeros(n, 1); % 相对位移 v zeros(n, 1); % 相对速度 a zeros(n, 1); % 相对加速度 p -m * ag(:); % 地震等效荷载 a(1) (p(1) - c * v(1) - k * u(1)) / m; for i 1 : n - 1 % 预测位移和速度 up u(i) dt * v(i) dt^2 * (0.5 - beta) * a(i); vp v(i) dt * (1 - gamma) * a(i); % 求解下一时刻相对加速度 a(i 1) (p(i 1) - c * vp - k * up) / ... (m gamma * dt * c beta * dt^2 * k); % 校正位移和速度 u(i 1) up beta * dt^2 * a(i 1); v(i 1) vp gamma * dt * a(i 1); end Sd max(abs(u)); % 相对位移谱 Sv max(abs(v)); % 相对速度谱 Sa max(abs(a ag(:))); % 绝对加速度谱 end代码里的预测-校正结构是 Newmark 格式的核心。先用上一时刻的位移、速度、加速度预测本时刻的 up 和 vp再把 up、vp 代进运动方程求出新的加速度最后用加速度校正位移和速度。这里的 p(i1)-m*ag(i1)单位与 ag 相同。返回的 Sa 写成 aag 而不是 a正是因为 2.1 节强调的绝对加速度定义。参数方面阻尼比 xi 按结构类型取钢筋混凝土结构通常取 0.05周期 T 越小 k 越大周期越大 k 越小。注意不要把单位混用如果 ag 用 cm/s^2则位移单位为 cm速度单位为 cm/s加速度单位为 cm/s^2后续与重力加速度 g981 对比时才一致。3. 反应谱计算参数怎么定阻尼比、周期布点与积分步长3.1 阻尼比按结构类型取不要统一用 0.05阻尼比是反应谱形状最敏感的参数。规范里混凝土结构取 0.05钢结构取 0.020.03隔震结构或耗能结构按附加阻尼取 0.10 甚至更高。阻尼比变大谱峰明显降低特别是靠近结构共振周期附近的峰值降幅最大。计算时把它做成函数参数而不是写死在代码里方便做参数敏感性分析也能直接对比不同阻尼假设对设计的差异。3.2 周期布点用 logspace不要用等间距数组反应谱横轴通常跨 0.02s 到 6s 甚至更长如果用 0.02:0.02:6 这种等间隔数组0.02s 到 1s 之间只有约 50 个点高频段谱型会非常粗糙峰值有可能被漏掉。常见做法是用对数均匀分布T_list logspace(log10(0.02), log10(6.0), 120);这条命令生成 0.02s 到 6s 之间按对数均匀分布的 120 个周期点。0.02s 到 0.1s 之间约有 20 个点1s 到 6s 之间也约有 20 个点都比等间距数组合理。点数太少曲线不够光滑点数太多计算时间线性增长120 到 240 个点是常用范围。结构设计复核时再在场地特征周期 Tg 附近加密几个点例如 Tg0.35s 时补 0.3、0.35、0.4 三个点能更好地读取平台段末端的峰值。这个加密习惯在比对规范设计谱时尤其有用。3.3 积分步长与最小可信周期Newmark-β 平均加速度法是无条件稳定的理论上可以直接使用地震记录的原始采样步长。但稳定性不等于精度当结构周期 T 小于 10 倍积分步长 dt 时高频段的数值色散会明显影响 Sv 和 Sa 峰值。比如 dt0.005s 的强震记录T0.02s 时一个周期只有 4 个积分步计算得到的相对速度谱可信度有限。所以程序里 T_min 默认取 0.02s如果记录采样率是 200Hz这已经是 4 个点作为工程估计够用但不建议把 T_min 压到 0.01s 以下。还有一个容易被忽略的输入处理环节是零线校正。实际强震记录在积分前需要先减去均值否则积分后的位移时程会产生明显漂移长周期段的 Sd 被显著高估。滤波则按需选用带通滤波但要注意高通截止频率取得太高会把长周期位移成分滤掉导致长周期 Sd 反而偏小。参数汇总如下参数建议取值说明阻尼比 xi混凝土 0.05钢结构 0.020.03谱峰高度随 xi 减小而增大周期范围0.026s科研场景可延伸到 10s周期点数120240 点对数均匀点数太少会错过谱峰积分步长直接用记录 dt0.005s无条件稳定但注意高频精度零线校正ag-mean(ag)必须做否则长周期段漂移提示如果改用 Duhamel 积分而不是 Newmark 隐式格式显式递推对步长有额外稳定性限制不要照搬这里的“直接用 dt 即可”结论。3.4 三个最容易被弄错的输入错误第一忽略记录的采样间隔直接用下标当时间轴导致 dt 写错一位整个谱曲线整体平移。第二把加速度单位与重力加速度单位混用在 cm/s^2 和 m/s^2 之间来回切换时不做换算结果 Sa 数量级差 100 倍几乎无法校验。第三反应谱计算必须用原始加速度记录逐点积分而不是先对加速度积分得到位移再对位移做傅里叶变换换算谱值后者得到的是谐波意义下的频响函数与瞬态地震反应的峰值没有直接对应关系。4. 完整流程从 csv 读入加速度记录到画出反应谱曲线4.1 用 readmatrix 读取 csv 并做零线校正data readmatrix(elcentro.csv); % 第一列时间(s)第二列加速度(gal) t data(:, 1); ag data(:, 2); ag(isnan(ag)) 0; % 空值置零避免污染后续积分 dt 0.005; % 期望采样间隔 dt_raw median(diff(t)); if abs(dt_raw - dt) / dt 0.01 % 与目标 dt 偏差超过 1% 时重新插值 tq 0 : dt : t(end); ag interp1(t, ag, tq, linear, 0); end ag ag(:) - mean(ag); % 零线校正readmatrix在较新版本的 Matlab 中可以直接读取带表头或不带表头的 csv。这里先检查时间序列的中值间隔与目标 dt 是否一致因为有些公开记录的采样率是 100Hz、200Hz 甚至 50Hz直接假设 0.005s 会出错。空值置零是保守处理如果空值较多建议回到原始数据重新检查而不是依赖置零。零线校正为什么必须做传感器即使标定良好也常有直流偏置直接带入运动方程会得到越来越大的速度漂移长周期段 Sd 完全失真。去均值是最廉价的处理。插值是自动判断而非强制因为强迫插值会改变原始高频成分。如果需要带通滤波fs 1 / dt; [b, a] butter(4, [0.1 25] / (fs / 2), bandpass); ag_f filtfilt(b, a, ag);需要信号处理工具箱。filtfilt是零相位滤波和filter的相位滞后有本质区别反应谱对相位敏感不要用filter替代。高通取 0.1Hz 会把长周期位移成分部分滤掉若你的关注点在最长周期的 Sd 值建议把高通降到 0.05Hz 或不做滤波。4.2 逐周期计算谱值并输出表格xi 0.05; T_list logspace(log10(0.02), log10(6), 120); nT length(T_list); Sd zeros(nT, 1); Sv zeros(nT, 1); Sa zeros(nT, 1); for j 1 : nT [Sd(j), Sv(j), Sa(j)] sdof_response_spectrum(ag_f, dt, xi, T_list(j)); end T_table table(T_list, Sd, Sv, Sa); writetable(T_table, response_spectrum.csv);外层循环在 120 个周期点、几千个采样点的情况下Matlab 计算时间在秒级不需要向量化优化。输出成 csv 表格便于后续比较不同阻尼比、不同记录的谱值也方便在 Excel 里做二次处理。如果对多个阻尼比循环计算记得每次循环重新分配结果数组不要在循环内用[Sd; new]动态拼接那会显著拖慢速度。4.3 用 loglog 绘制三条谱曲线并标注单位反应谱横纵轴跨度都大线性坐标下低频段和高频段无法同时看清。使用loglogfigure(Color, w); loglog(T_list, Sd, b-, LineWidth, 1.2); hold on; loglog(T_list, Sv, r-, LineWidth, 1.2); loglog(T_list, Sa, k-, LineWidth, 1.2); hold off; legend(相对位移谱 Sd (cm), ... 相对速度谱 Sv (cm/s), ... 绝对加速度谱 Sa (cm/s^2), Location, northwest); xlabel(周期 T (s)); ylabel(反应谱值); grid on;三条线的量纲不同所以 legend 里必须带单位。图形采用双对数坐标后短周期段的细节能拉开长周期段的衰减趋势也更接近工程习惯。如果之后要把绝对加速度谱和设计谱叠加单位需要统一换算成无量纲地震影响系数 α这个放在第 5 章处理。4.4 用一个正弦输入对拍程序正确性最好的自查方式是用一条已知解析解的正弦激励对拍。取周期 T2s 的单自由度体系输入 0.5Hz、幅值 0.1g 的正弦波阻尼比取 0.01。在共振稳态下位移和绝对加速度的放大倍数近似为 1/(2ξ)50因此 Sa 应接近 0.1×505g即 490 cm/s^2。fs 200; dt 1 / fs; t (0 : dt : 40); ag_test 0.1 * 981 * sin(2 * pi * 0.5 * t); % 0.5Hz, 0.1g ag_test(1:round(5/dt)) 0; % 前 5s 置零等待稳态建立 [Sd, Sv, Sa] sdof_response_spectrum(ag_test, dt, 0.01, 2.0);预期 Sa 约在 490 cm/s^2 附近误差超过 5% 时优先检查阻尼比和周期是否传错再看输入波形幅值是否包含了 981 这个重力加速度换算。注意前 5s 置零是为了让体系从静止开始过渡到稳态避免瞬态响应对峰值的干扰如果不加这段置零Sa 会包含启动瞬态分量与稳态放大倍数的解析解对不上。5. 用绝对加速度谱快速估算底部剪力并完成三道自检5.1 从 Sa 换算地震影响系数 α工程设计中关心的是绝对加速度谱 Sa 与重力加速度 g 的比值即地震影响系数 αSa/g。从已计算的离散谱点读取任意周期的 α 时建议在对数坐标系插值因为谱曲线在双对数坐标下更接近分段线性g 981; % cm/s^2 alpha Sa / g; T1 1.0; % 目标周期须落于 T_list 范围内 log_alpha interp1(log10(T_list), log10(alpha), log10(T1)); alpha1 10^log_alpha;单自由度结构底部剪力可写成 Vb α1 · m · g其中 m 为结构质量。这里的 α1 已经是考虑了阻尼比后的反应谱值可以直接用于弹性阶段的粗略估算。注意 T1 若超出 T_list 的范围interp1会返回 NaN因此不要读取范围之外的周期。对于长周期结构要么扩展 T_list 的上限要么用外插函数单独处理。5.2 三道自检高频端、长周期端、伪谱一致性程序跑完后别急着交付结果做三次快速检查就能发现八成以上的错误。第一高频端 Sa(T_min) 应该接近输入的峰值加速度 PGA。刚体极限下结构相对位移为零绝对加速度等于地面加速度因此最小周期处的 Sa 应与 PGA 同量级。如果滤波后的 PGA 与原始 PGA 差异明显就用滤波后的峰值比对。第二长周期端 Sd(T_max) 应与地面位移时程峰值的数量级一致比如积分地面加速度得到的位移峰值如果只有几厘米Sd 却报了上百厘米说明零线校正没做或高通滤波参数不对。第三伪速度谱 ω·Sd 与真实 Sv 在两个数量级范围内应基本重合PSv 2 * pi ./ T_list .* Sd; loglog(T_list, Sv, T_list, PSv, --); legend(真实 Sv,伪速度谱 \omega·Sd, Location, best);小阻尼ξ0.05下伪速度谱与真实 Sv 的偏差通常在 5%10% 以内且真实 Sv 在共振峰附近略低。若两条线在中短周期段严重分离优先怀疑阻尼比是否写成了 0.5或者积分步长与采样率不一致。这三道自检全部通过后输出的相对位移谱、相对速度谱和绝对加速度谱才具备可信度也才适合用作抗震分析结论的支撑材料。本文还有配套的精品资源点击获取 用Matlab实现地震反应谱计算：相对位移、速度与绝对加速度谱
简介针对地震工程分析需求这份MATLAB程序基于《抗震工程学》理论编写面向科研人员与在校学生可快速由地震波原始数据生成相对速度谱、相对位移谱和绝对加速度谱。资源包内共2个文件m源码文件用于执行谱计算与绘图txt文件存放支持运行的地震波数据整体压缩包仅11KB轻量易用。目前已有463人学习下载程序调试完备、精度可靠能够满足地震分析、课程作业及入门实践等场景。完成下载后即可直接加载数据运行通过结果直观对比三类反应谱曲线帮助理解结构地震响应特性也为后续研究提供了可扩展的代码基础和真实数据支撑。1. 同一张反应谱图为什么要分相对位移、相对速度和绝对加速度三条线拿到一条采样间隔 0.005s 的强震加速度记录时结构工程师通常不先问 PGA 是多少而是问周期 1s 左右的框架这条波能把它推到多大力。直接看波形回答不了——同一条波对周期 0.2s 和 2s 的单自由度体系的放大效应可能差一个数量级。反应谱就是为这个场景设计的对一组周期逐次求解单自由度运动方程提取相对位移、相对速度、绝对加速度的时程峰值得到相对位移谱 Sd、相对速度谱 Sv、绝对加速度谱 Sa 三条谱曲线。这里给出的 Matlab 程序覆盖从运动方程、Newmark 积分、参数标定到画图验证的完整链路适合做抗震分析、结构设计复核或课程作业的读者直接参照落地。2. 从运动方程到三条谱曲线反应谱在算什么2.1 单自由度体系的相对运动方程要理解相对位移谱跟绝对加速度谱的差别关键是把“相对”和“绝对”两个坐标系分清。取结构底部跟着地面一起平移结构质量 m 相对地面的位移记作 u(t)地面自身位移记作 ug(t)那么质量块的绝对位移是 u(t)ug(t)。以相对位移为未知量单自由度体系的运动方程为$$m\ddot{u}(t)c\dot{u}(t)ku(t)-m\ddot{u}_g(t)$$右侧的 $\ddot{u}_g(t)$ 就是地震记录里的加速度时程可以直接从 csv 或文本文件读进来。它不是传统意义上的外荷载而是地面运动引起的等效惯性力。这里的阻尼 c 和刚度 k 都不是随意给定的刚度 k 对应结构周期$km(2\pi/T)^2$阻尼 c 对应阻尼比$c2m\xi\omega$。实际程序里通常把质量 m 归一化成 1只保留阻尼比 $\xi$ 和圆频率 $\omega$ 两个自由度。绝对加速度在这里不是 $\ddot{u}(t)$而是 $\ddot{u}(t)\ddot{u}_g(t)$也就是质点相对地面的加速度加上地面本身的加速度。很多初学者以为反应谱里的“绝对加速度谱”指的是相对加速度的峰值这会导致后期跟规范里的地震影响系数对不上需要格外留意。2.2 相对位移谱、相对速度谱与绝对加速度谱的定义程序要输出的三条谱在文献里有固定的英文缩写和计算方式Sd 是相对位移时程绝对值的最大值Sv 是相对速度时程绝对值的最大值Sa 是绝对加速度时程绝对值的最大值。对同一条地震波不同的周期 T 会得到不同的 Sd、Sv、Sa 值把所有这些值连成曲线就构成反应谱。工程上还有一个经常出现的计算变体就是伪速度谱和伪加速度谱分别定义为 $\omega S_d$ 和 $\omega^2 S_d$。教科书和规范里的速度反应谱多数是伪速度谱而不是真实相对速度谱。两者在小阻尼下数值非常接近但定义不同。本程序的默认输出是严格定义值$S_d\max_t|u(t)|$$S_v\max_t|\dot u(t)|$$S_a\max_t|\ddot u(t)\ddot u_g(t)|$同时把伪谱作为校验手段保留在验证环节。2.3 Newmark-β 平均加速度法的 Matlab 实现时域积分使用 Newmark-β 法里最常用的平均加速度格式$\gamma1/2$、$\beta1/4$。它无条件稳定可以使用原始采样步长而不必为了稳定性加密积分步长。下面是单自由度体系反应谱计算函数输入一条加速度时程、采样间隔、阻尼比和周期输出三个谱值。function [Sd, Sv, Sa] sdof_response_spectrum(ag, dt, xi, T) % ag : 地面加速度时程单位任意cm/s^2 或 m/s^2 % dt : 采样间隔单位 s % xi : 阻尼比例如 0.05 % T : 单自由度体系周期单位 s % Sd Sv Sa 分别为相对位移、相对速度、绝对加速度峰值 omega 2 * pi / T; % 圆频率 m 1; % 质量归一化 c 2 * xi * omega * m; % 阻尼系数 k omega^2 * m; % 刚度系数 % Newmark 平均加速度常数 gamma 1 / 2; beta 1 / 4; n length(ag); u zeros(n, 1); % 相对位移 v zeros(n, 1); % 相对速度 a zeros(n, 1); % 相对加速度 p -m * ag(:); % 地震等效荷载 a(1) (p(1) - c * v(1) - k * u(1)) / m; for i 1 : n - 1 % 预测位移和速度 up u(i) dt * v(i) dt^2 * (0.5 - beta) * a(i); vp v(i) dt * (1 - gamma) * a(i); % 求解下一时刻相对加速度 a(i 1) (p(i 1) - c * vp - k * up) / ... (m gamma * dt * c beta * dt^2 * k); % 校正位移和速度 u(i 1) up beta * dt^2 * a(i 1); v(i 1) vp gamma * dt * a(i 1); end Sd max(abs(u)); % 相对位移谱 Sv max(abs(v)); % 相对速度谱 Sa max(abs(a ag(:))); % 绝对加速度谱 end代码里的预测-校正结构是 Newmark 格式的核心。先用上一时刻的位移、速度、加速度预测本时刻的 up 和 vp再把 up、vp 代进运动方程求出新的加速度最后用加速度校正位移和速度。这里的 p(i1)-m*ag(i1)单位与 ag 相同。返回的 Sa 写成 aag 而不是 a正是因为 2.1 节强调的绝对加速度定义。参数方面阻尼比 xi 按结构类型取钢筋混凝土结构通常取 0.05周期 T 越小 k 越大周期越大 k 越小。注意不要把单位混用如果 ag 用 cm/s^2则位移单位为 cm速度单位为 cm/s加速度单位为 cm/s^2后续与重力加速度 g981 对比时才一致。3. 反应谱计算参数怎么定阻尼比、周期布点与积分步长3.1 阻尼比按结构类型取不要统一用 0.05阻尼比是反应谱形状最敏感的参数。规范里混凝土结构取 0.05钢结构取 0.020.03隔震结构或耗能结构按附加阻尼取 0.10 甚至更高。阻尼比变大谱峰明显降低特别是靠近结构共振周期附近的峰值降幅最大。计算时把它做成函数参数而不是写死在代码里方便做参数敏感性分析也能直接对比不同阻尼假设对设计的差异。3.2 周期布点用 logspace不要用等间距数组反应谱横轴通常跨 0.02s 到 6s 甚至更长如果用 0.02:0.02:6 这种等间隔数组0.02s 到 1s 之间只有约 50 个点高频段谱型会非常粗糙峰值有可能被漏掉。常见做法是用对数均匀分布T_list logspace(log10(0.02), log10(6.0), 120);这条命令生成 0.02s 到 6s 之间按对数均匀分布的 120 个周期点。0.02s 到 0.1s 之间约有 20 个点1s 到 6s 之间也约有 20 个点都比等间距数组合理。点数太少曲线不够光滑点数太多计算时间线性增长120 到 240 个点是常用范围。结构设计复核时再在场地特征周期 Tg 附近加密几个点例如 Tg0.35s 时补 0.3、0.35、0.4 三个点能更好地读取平台段末端的峰值。这个加密习惯在比对规范设计谱时尤其有用。3.3 积分步长与最小可信周期Newmark-β 平均加速度法是无条件稳定的理论上可以直接使用地震记录的原始采样步长。但稳定性不等于精度当结构周期 T 小于 10 倍积分步长 dt 时高频段的数值色散会明显影响 Sv 和 Sa 峰值。比如 dt0.005s 的强震记录T0.02s 时一个周期只有 4 个积分步计算得到的相对速度谱可信度有限。所以程序里 T_min 默认取 0.02s如果记录采样率是 200Hz这已经是 4 个点作为工程估计够用但不建议把 T_min 压到 0.01s 以下。还有一个容易被忽略的输入处理环节是零线校正。实际强震记录在积分前需要先减去均值否则积分后的位移时程会产生明显漂移长周期段的 Sd 被显著高估。滤波则按需选用带通滤波但要注意高通截止频率取得太高会把长周期位移成分滤掉导致长周期 Sd 反而偏小。参数汇总如下参数建议取值说明阻尼比 xi混凝土 0.05钢结构 0.020.03谱峰高度随 xi 减小而增大周期范围0.026s科研场景可延伸到 10s周期点数120240 点对数均匀点数太少会错过谱峰积分步长直接用记录 dt0.005s无条件稳定但注意高频精度零线校正ag-mean(ag)必须做否则长周期段漂移提示如果改用 Duhamel 积分而不是 Newmark 隐式格式显式递推对步长有额外稳定性限制不要照搬这里的“直接用 dt 即可”结论。3.4 三个最容易被弄错的输入错误第一忽略记录的采样间隔直接用下标当时间轴导致 dt 写错一位整个谱曲线整体平移。第二把加速度单位与重力加速度单位混用在 cm/s^2 和 m/s^2 之间来回切换时不做换算结果 Sa 数量级差 100 倍几乎无法校验。第三反应谱计算必须用原始加速度记录逐点积分而不是先对加速度积分得到位移再对位移做傅里叶变换换算谱值后者得到的是谐波意义下的频响函数与瞬态地震反应的峰值没有直接对应关系。4. 完整流程从 csv 读入加速度记录到画出反应谱曲线4.1 用 readmatrix 读取 csv 并做零线校正data readmatrix(elcentro.csv); % 第一列时间(s)第二列加速度(gal) t data(:, 1); ag data(:, 2); ag(isnan(ag)) 0; % 空值置零避免污染后续积分 dt 0.005; % 期望采样间隔 dt_raw median(diff(t)); if abs(dt_raw - dt) / dt 0.01 % 与目标 dt 偏差超过 1% 时重新插值 tq 0 : dt : t(end); ag interp1(t, ag, tq, linear, 0); end ag ag(:) - mean(ag); % 零线校正readmatrix在较新版本的 Matlab 中可以直接读取带表头或不带表头的 csv。这里先检查时间序列的中值间隔与目标 dt 是否一致因为有些公开记录的采样率是 100Hz、200Hz 甚至 50Hz直接假设 0.005s 会出错。空值置零是保守处理如果空值较多建议回到原始数据重新检查而不是依赖置零。零线校正为什么必须做传感器即使标定良好也常有直流偏置直接带入运动方程会得到越来越大的速度漂移长周期段 Sd 完全失真。去均值是最廉价的处理。插值是自动判断而非强制因为强迫插值会改变原始高频成分。如果需要带通滤波fs 1 / dt; [b, a] butter(4, [0.1 25] / (fs / 2), bandpass); ag_f filtfilt(b, a, ag);需要信号处理工具箱。filtfilt是零相位滤波和filter的相位滞后有本质区别反应谱对相位敏感不要用filter替代。高通取 0.1Hz 会把长周期位移成分部分滤掉若你的关注点在最长周期的 Sd 值建议把高通降到 0.05Hz 或不做滤波。4.2 逐周期计算谱值并输出表格xi 0.05; T_list logspace(log10(0.02), log10(6), 120); nT length(T_list); Sd zeros(nT, 1); Sv zeros(nT, 1); Sa zeros(nT, 1); for j 1 : nT [Sd(j), Sv(j), Sa(j)] sdof_response_spectrum(ag_f, dt, xi, T_list(j)); end T_table table(T_list, Sd, Sv, Sa); writetable(T_table, response_spectrum.csv);外层循环在 120 个周期点、几千个采样点的情况下Matlab 计算时间在秒级不需要向量化优化。输出成 csv 表格便于后续比较不同阻尼比、不同记录的谱值也方便在 Excel 里做二次处理。如果对多个阻尼比循环计算记得每次循环重新分配结果数组不要在循环内用[Sd; new]动态拼接那会显著拖慢速度。4.3 用 loglog 绘制三条谱曲线并标注单位反应谱横纵轴跨度都大线性坐标下低频段和高频段无法同时看清。使用loglogfigure(Color, w); loglog(T_list, Sd, b-, LineWidth, 1.2); hold on; loglog(T_list, Sv, r-, LineWidth, 1.2); loglog(T_list, Sa, k-, LineWidth, 1.2); hold off; legend(相对位移谱 Sd (cm), ... 相对速度谱 Sv (cm/s), ... 绝对加速度谱 Sa (cm/s^2), Location, northwest); xlabel(周期 T (s)); ylabel(反应谱值); grid on;三条线的量纲不同所以 legend 里必须带单位。图形采用双对数坐标后短周期段的细节能拉开长周期段的衰减趋势也更接近工程习惯。如果之后要把绝对加速度谱和设计谱叠加单位需要统一换算成无量纲地震影响系数 α这个放在第 5 章处理。4.4 用一个正弦输入对拍程序正确性最好的自查方式是用一条已知解析解的正弦激励对拍。取周期 T2s 的单自由度体系输入 0.5Hz、幅值 0.1g 的正弦波阻尼比取 0.01。在共振稳态下位移和绝对加速度的放大倍数近似为 1/(2ξ)50因此 Sa 应接近 0.1×505g即 490 cm/s^2。fs 200; dt 1 / fs; t (0 : dt : 40); ag_test 0.1 * 981 * sin(2 * pi * 0.5 * t); % 0.5Hz, 0.1g ag_test(1:round(5/dt)) 0; % 前 5s 置零等待稳态建立 [Sd, Sv, Sa] sdof_response_spectrum(ag_test, dt, 0.01, 2.0);预期 Sa 约在 490 cm/s^2 附近误差超过 5% 时优先检查阻尼比和周期是否传错再看输入波形幅值是否包含了 981 这个重力加速度换算。注意前 5s 置零是为了让体系从静止开始过渡到稳态避免瞬态响应对峰值的干扰如果不加这段置零Sa 会包含启动瞬态分量与稳态放大倍数的解析解对不上。5. 用绝对加速度谱快速估算底部剪力并完成三道自检5.1 从 Sa 换算地震影响系数 α工程设计中关心的是绝对加速度谱 Sa 与重力加速度 g 的比值即地震影响系数 αSa/g。从已计算的离散谱点读取任意周期的 α 时建议在对数坐标系插值因为谱曲线在双对数坐标下更接近分段线性g 981; % cm/s^2 alpha Sa / g; T1 1.0; % 目标周期须落于 T_list 范围内 log_alpha interp1(log10(T_list), log10(alpha), log10(T1)); alpha1 10^log_alpha;单自由度结构底部剪力可写成 Vb α1 · m · g其中 m 为结构质量。这里的 α1 已经是考虑了阻尼比后的反应谱值可以直接用于弹性阶段的粗略估算。注意 T1 若超出 T_list 的范围interp1会返回 NaN因此不要读取范围之外的周期。对于长周期结构要么扩展 T_list 的上限要么用外插函数单独处理。5.2 三道自检高频端、长周期端、伪谱一致性程序跑完后别急着交付结果做三次快速检查就能发现八成以上的错误。第一高频端 Sa(T_min) 应该接近输入的峰值加速度 PGA。刚体极限下结构相对位移为零绝对加速度等于地面加速度因此最小周期处的 Sa 应与 PGA 同量级。如果滤波后的 PGA 与原始 PGA 差异明显就用滤波后的峰值比对。第二长周期端 Sd(T_max) 应与地面位移时程峰值的数量级一致比如积分地面加速度得到的位移峰值如果只有几厘米Sd 却报了上百厘米说明零线校正没做或高通滤波参数不对。第三伪速度谱 ω·Sd 与真实 Sv 在两个数量级范围内应基本重合PSv 2 * pi ./ T_list .* Sd; loglog(T_list, Sv, T_list, PSv, --); legend(真实 Sv,伪速度谱 \omega·Sd, Location, best);小阻尼ξ0.05下伪速度谱与真实 Sv 的偏差通常在 5%10% 以内且真实 Sv 在共振峰附近略低。若两条线在中短周期段严重分离优先怀疑阻尼比是否写成了 0.5或者积分步长与采样率不一致。这三道自检全部通过后输出的相对位移谱、相对速度谱和绝对加速度谱才具备可信度也才适合用作抗震分析结论的支撑材料。本文还有配套的精品资源点击获取