简介一份面向结构工程、地震工程初学者与MATLAB使用者的完整数据处理资源包围绕从PEER地震研究中心获取实测地震波数据基于杜哈梅积分将加速度时程转换为位移响应并通过FFT频域变换分析频谱特征系统讲解了振动信号处理与结构动力响应计算流程。资源共913个文件压缩包约20.51MB主体为三类各302个的at2、vt2、dt2格式地震记录文件另有3个MATLAB脚本Duhamel.m、FFT.m等、2份PPT教学课件、1个xls及1个csv数据表便于对照代码理解每个计算环节。已有3121人学习使用适合正在做振动信号分析作业或地震波处理相关课题的学生参考。资源内附PPT解释了理论公式与实现步骤通过NORTHRIDGE地震实例演示从PEER原始记录到位移曲线的完整链路既能掌握杜哈梅积分数值实现也能学会用MATLAB进行频域特征提取是一份理论与实践结合紧密的辅助材料。1. 地震波进了Matlab第一步不是积分而是先问数据从哪来做结构抗震分析的人多半都卡过这一步从PEER NGA数据库筛出一条符合场地和震级的地震波下载下来一看是AT2格式直接扔进Matlab里plot发现全乱了。更常见的情况是加速度记录进来了但你要的位移时程怎么算都跟参考曲线对不上频域结果画出来还一头是尖刺一头是噪声。这个标题把Matlab、PEER、duhamel积分和FFT四个词拼在一起本质上是在讲一条完整的地震波处理链路下载格式解析、加速度基线校正、时域杜哈梅积分求位移、再用快速傅里叶变换从频域交叉验证。这套流程在本科毕业设计、研究生课题和实际工程弹性时程分析里都能直接用是结构动力学课里“拿真实数据验证理论公式”的典型场景。本文按我平时处理这类数据的顺序展开先把PEER下载的AT2文件解析成标准列向量再写一个不依赖内置函数的duhamel积分器因为教材公式和实际离散采样之间差着好几层坑然后用FFT做频域求解做对照最后把两条位移曲线画叠在一块验证一致性。整个过程会给出可以直接改改就用的Matlab代码并解释每一步的关键参数。2. 从PEER下载到MatlabAT2格式解析与单位陷阱2.1 AT2文件长什么样为什么Matlab的load读不了PEER NGA-West2数据库是结构工程领域最常用的强震记录来源。下载单条记录时会得到一个逐行文本文件后缀是AT2内部结构大致分三层文件头若干行、加速度数据行、末尾可能带空行。文件头里藏着关键信息比如时间步长DT、点数NPTS、地震事件名称、台站名称、震级和震中距。真正让Matlab的load函数无能为力的是数据段采用了固定列宽的Fortran格式输出每行通常是5个浮点数每个字段占15个字符左右字段间没有逗号也没有制表符。如果你直接读文件再split会发现数字粘在一起或者小数点位错乱。这就是为什么必须用textread或fscanf按固定格式解析。读取之前还要看一眼文件头里的UNITS行PEER默认给出的加速度单位是g1g约等于9.80665 m/s²如果后面要做duhamel积分得出以米为单位的位移就得先做单位换算。很多人在这一步漏了换算导致最终位移曲线数值大了快十倍。2.2 用fgetl和sscanf把AT2变成Matlab列向量最稳妥的做法是用fgetl逐行读文件头跳过不关心的行然后从数据段起始行开始用sscanf按浮点数批量解析。下面这段代码直接处理PEER标准AT2文件并返回时间轴、加速度m/s²和时间步长。function [t, acc, dt] read_peer_at2(filename) % 读取PEER NGA-West2的AT2格式加速度记录 % 返回 % t - 时间向量 (s) % acc - 加速度向量 (m/s^2) % dt - 采样时间间隔 (s) fid fopen(filename, r); if fid -1 error(无法打开文件: %s, filename); end % 读前几行定位数据起始 % 注意不同版本的PEER文件头行数可能不同 % 所以用while循环找数据行而不是硬编码行号 data_started false; header_line ; np 0; dt 0.01; % 默认值后面会被文件头中的真实值覆盖 while ~feof(fid) line fgetl(fid); % 文件头中的关键标记点数NPTS和时间步长DT if contains(line, NPTS) vals sscanf(line, %*s %*s %d %*s %*s %f); if length(vals) 2 np vals(1); dt vals(2); end end if contains(line, ACCELERATION) data_started true; % 跳过“ACCELERATION”这一行本身 continue; end if data_started header_line line; break; end end % 现在header_line是第一个数据行但后面还有更多行 % 把当前行和剩余行全部读进来 data_str header_line; while ~feof(fid) data_str [data_str, , fgetl(fid)]; %#okAGROW end fclose(fid); % 用sscanf按浮点数解析 raw sscanf(data_str, %f); % 取前np个点防止末尾空行混入脏数据 acc raw(1:np); % PEER默认单位为g换算为m/s^2 acc acc * 9.80665; t (0:np-1) * dt; end这段代码里最关键的是sscanf(data_str, %f)它不关心字段占多少列只按连续数值逐块读取天然适配AT2的固定列宽布局。另一点是单位换算放在了解析之后、积分之前避免在duhamel积分结果里引入量纲错误。如果文件头里没有ACCELERATION标记说明文件格式不是标准PEER布局此时应当打印前20行人工确认数据起始位置。2.3 特别提示先画时程图再往下走拿到加速度后第一件事不是立刻积分而是plot出来看形态。真实地震记录的首尾会有偏移中间有几段明显的高频振动是正常的但如果你看到一条斜率恒定的直线说明基线没归零后面duhamel积分出来的位移会飞掉。此外还要确认时间轴总长度PEER记录通常截断在有效振动段之后有些记录只有15秒有些长到60秒以上这直接影响后续积分步数和FFT频率分辨率。遇到加速度首尾不为零的基线漂移时常见做法是对整条记录减去均值更讲究的会用高通滤波去掉低频漂移。但滤波会改变波形本身在弹性时程分析里一般只做均值修正就够了。我的习惯是保存原始加速度向量再单独保存修正后的向量后面做FFT对比时两种都画避免调试时搞不清曲线差异是算法问题还是数据预处理问题。3. 杜哈梅积分离散化从教材公式到能跑的代码3.1 为什么不能直接对加速度积分两次很多初学者试图对加速度做两次cumtrapz得到位移这在数值数学上没错但在地震工程里是致命的。原因有两个真实地震加速度记录含有低频噪声和仪器漂移双重积分会把低频误差放大成抛物线趋势位移曲线末端可能离谱到几十米结构响应是动力放大过程简单积分完全忽略体系的刚度、阻尼和自振周期特征得到的是地面运动位移不是结构相对位移响应。杜哈梅积分的意义在于它描述了单自由度体系在地面加速度作用下的动力响应。对线性单自由度体系相对位移响应表达式是一个卷积积分也就是把地震加速度记录与体系的脉冲响应函数做卷积。教材里那个带阻尼的积分表达式看起来很干净但真实地震记录是离散采样的必须做数值离散化处理这里就涉及积分步长和阻尼项的处理直接套教材公式会跑出振荡发散的结果。3.2 离散化的三种常见递推方案工程中最常用的离散化递推格式有Newmark-β法、中心差分法和分段精确积分法Duhamel积分的解析离散形式。中心差分法是显式格式计算量小但稳定条件要求时间步长足够小PEER记录采样率通常是100Hz或200Hz即dt为0.01秒或0.005秒一般能满足要求。Newmark-β法是隐式格式无条件稳定适合自振周期较长的结构。但如果是“做duhamel积分”为主题的代码我更倾向于直接用分段精确积分法对每个时间步内假设加速度线性变化可以推导出递推公式精度高于中心差分也不需要像Newmark那样选择β和γ参数。下面是阻尼单自由度体系的分段精确积分递推实现function [u, v, a] duhamel_response(t, acc, Tn, zeta) % 分段精确积分法求解单自由度体系地震响应 % 输入 % t - 时间向量 (s) % acc - 地面运动加速度 (m/s^2) % Tn - 结构自振周期 (s) % zeta - 阻尼比 (0.01~0.05, 一般取0.05) % 输出 % u - 相对位移响应 (m) % v - 相对速度响应 (m/s) % a - 相对加速度响应 (m/s^2) n length(acc); u zeros(n, 1); v zeros(n, 1); a zeros(n, 1); dt t(2) - t(1); wn 2 * pi / Tn; wd wn * sqrt(1 - zeta^2); % 指数衰减项系数 e1 exp(-zeta * wn * dt); c cos(wd * dt); s sin(wd * dt); % 递推系数矩阵由状态空间解析解导出 A11 e1 * (c zeta * wn / wd * s); A12 e1 * (s / wd); A21 -e1 * (wn^2 / wd * s); A22 e1 * (c - zeta * wn / wd * s); % 激励项系数加速度在时段内线性变化 B11 e1 * ((zeta * wn / (wn^2) * c (1/wd - zeta*wn/(wn^2*wd)) * s) / 1 - 1/wn^2); % 实际工程里一般把地面加速度当成片段常数下面用简化系数 h dt; % 对每个时间步做递推 for i 1:n-1 % 地面加速度在步内取平均值 p_i acc(i); p_ip1 acc(i1); % 线性变化激励的解析间断解系数直接展开 coeff (p_ip1 - p_i) / dt; term1 A11 * u(i) A12 * v(i); term2 (p_i / wn^2) * (1 - A11 - zeta*wn*A12); term3 (coeff / wn^2) * (dt - A12 - zeta*wn*(dt - A11)); u(i1) term1 term2 term3; v(i1) A21 * u(i) A22 * v(i) ... (p_i / wn^2) * (zeta*wn*(1 - A11 - zeta*wn*A12) - wn^2*A12) ... (coeff / wn^2) * (1 - A21 - zeta*wn*(1 - A11) - wn^2*A12); a(i1) -2 * zeta * wn * v(i1) - wn^2 * u(i1) - acc(i1); end end说明一下这段代码的设计逻辑状态空间离散化得到的是解析递推解步内假设加速度线性变化这比直接数值积分精度高很多。a(i1)的计算用的是运动方程也就是相对加速度等于负的阻尼力减去弹性力再减去地面加速度。如果你想校验递推系数的正确性可以把zeta设为0、Tn设得很大此时结构退化为刚体位移应接近地面位移的负值用这个极限场景做验证非常有效。3.3 四个必调的参数杜哈梅积分器运行前要确认四个参数结构自振周期Tn、阻尼比zeta、采样时间间隔dt和积分总时长。前三者直接决定递推矩阵的数值特性dt不能大于结构自振周期的十分之一否则即便递推格式是稳的计算结果也会失真。阻尼比做弹性分析时取0.05是规范惯例但如果你想观察不同阻尼对位移峰值的影响可以写成循环扫描0.02到0.1。总时长方面如果地震记录只有15秒而结构自振周期是2秒那么至少得保证记录后段加速度衰减到接近零否则截断效应会让尾部位移出现不应有的振动。一个简单的检查方式是看加速度时程最后5秒钟的均方根值是否小于峰值的5%。如果没达到你需要在记录后面补零让结构有时间完成自由衰减振动。补零的操作在duhamel积分里是合法的因为那代表地面运动已经停止。3.4 结果合理性判断跑完积分后先画出位移时程肉眼检查形态初始段接近零中间段随强震到达而增大末尾段做衰减振荡。如果末端出现线性发散大概率是基线漂移没有修正回到第2步重新做均值归零。如果曲线整体幅值极小比如最大位移只有几毫米检查是不是单位换算出错加速度还停留在g。如果曲线呈现明显的高频毛刺检查dt和结构周期的比值是否太小。还有一种常见现象是位移曲线峰值出现在强震到达瞬间附近但略微滞后这是阻尼体系对激励的相位滞后是正常物理现象不要当作bug去改代码。4. FFT频域求解位移曲线用卷积定理做交叉验证4.1 为什么有了时域结果还要做FFT时域duhamel积分直接给出位移曲线那FFT频域求解看起来像是绕远路但它在验证算法正确性上价值很高。线性体系在频域里的传递函数是明确的代数表达式地面加速度做FFT得到频谱乘以传递函数得到位移频谱再做IFFT返回时域整个过程里没有任何递推累计误差。两条曲线如果对得上说明时域递推系数矩阵没有写错积分公式无误。这个交叉验证框架也是工程上常用的做法两个独立算法结果互相印证比单独一次计算更可信。更进一步FFT求解过程能同时输出位移频谱方便看到结构自振频率附近的放大效应这是时域曲线里不容易直接读出来的信息。4.2 构造传递函数单自由度体系频响函数 H(ω)线性单自由度体系在频域中满足代数方程相对位移对地面加速度的频响函数为H(ω) -1 / (ωn² - ω² 2 * i * ζ * ωn * ω)其中ωn是结构自振圆频率ω是激励频率。把这个解析表达式用Matlab代码写出来配合fft函数就能完成整个频域求解链路。要注意的是Matlab的fft输出是单边还是双边的问题不影响计算正确性只要你保证乘法在频域逐点对应、最后用ifft时除以点数幅值关系自然正确。下面这段代码是频域求解的核心function u_freq fft_response(t, acc, Tn, zeta) % 频域法求解单自由度体系位移响应 % 输入输出含义与duhamel_response保持一致 n length(acc); dt t(2) - t(1); Fs 1 / dt; % 加速度信号做FFT ACC fft(acc); % 构建频率轴注意fft频率分布是0到Fs然后负频率 f (0:n-1) * Fs / n; w 2 * pi * f; % 结构参数 wn 2 * pi / Tn; % 频响函数 H(w) -1 / (wn^2 - w^2 2i*zeta*wn*w) % 注意这里的负号源自运动方程中地面加速度项的符号 H -1 ./ (wn^2 - w.^2 2i * zeta * wn * w); % 频域相乘得到位移频谱 U ACC .* H; % IFFT返回时域 u_freq real(ifft(U)); % 因为是双边谱IFFT后幅值直接是真实响应不需要额外缩放 end这段代码里最容易出错的点在频响函数分母的符号。如果你把阻尼项写成了2i*zeta*wn*w前面是负号在时域里会变成负阻尼曲线会越振越大。另一个细节是w包含了0频率分量也就是直流项对于真实地震波而言直流分量极小但在公式计算中0频率处分母为wn²不会导致除零错误。若结构自振周期极大导致wn接近0才会出现除法溢出那时候要加一个eps保护项。4.3 两个算法结果对不齐时看什么使用时域和频域各算一遍后把两条曲线画在同一张图上设置不同的线型便于区分。正常情况下两条曲线几乎完全重合峰值相对误差在1%以内。如果对不齐优先检查三个位置起始段、峰值段和结尾段。起始段对不齐说明初始条件设置不一致频域法默认体系从静止开始时域法递推初始条件也是零通常一致峰值段对不齐多半是阻尼输入不一致确认两个函数里zeta传入的是同一个数值结尾段对不齐往往是记录截断效应频域FFT隐含了周期性延拓假设而真实记录末尾加速度不为零时这种周期延拓会引入泄漏误差。解决FFT尾端泄漏的办法是给加速度记录加窗函数比如Hann窗或指数窗。但要注意加窗会改变时程形状导致峰值被压低所以加窗只用于频谱观察不用于最终位移计算。正确做法是把原始记录末尾补零到合适的长度再做FFT补零不改变信号内容只提高频率采样密度配合整周期截取能大幅减少泄漏。补零长度一般选择超过原始长度的下一个2的幂次比如原始1500点补到2048点。4.4 顺手把加速度频谱也画出来在验证位移曲线的同时我还习惯把加速度的傅里叶幅值谱画出来横轴用Hz纵轴用log刻度。这张图能直观看到地震波能量集中的频段合理选择结构自振周期时很有用。若把位移响应频谱也画上去会看到位移谱峰值出现在结构自振频率附近这正是共振现象的体现。这个双谱图可以作为一份报告里的中间插图用来支撑“为什么要选择这个结构周期”的结论。代码实现很短plot(f(1:n/2), abs(ACC(1:n/2)))取前半段是因为后半段是负频率共轭部分对实信号没有额外信息。横轴单位是Hz如果你读文件时dt是0.01秒那么Fs100HzNyquist频率是50Hz超过50Hz的频段根本没有有效信息别把横轴画到100Hz再讨论物理意义。5. 把两条曲线叠画在一张图上的验证脚本5.1 主脚本如何组织这些函数前面四个函数各自独立但一个完整的处理流程需要主脚本来调度。我的组织方式是用一个m文件做入口依次执行读取AT2、画加速度时程、时域积分、频域求解、画位移对比图、画频谱图。这样改文件名和参数时只动主脚本函数内部逻辑不用碰。主脚本里定义结构体存储地震信息方便最后统一输出地震动峰值加速度PGA和最大位移响应。% 主脚本: run_seismic_analysis.m % 使用前请修改为你的AT2文件实际路径 filename RSN1234_EQNAME_STATION.AT2; % 1. 读取PEER记录 [t, acc, dt] read_peer_at2(filename); % 2. 基线修正去均值 acc_detrend acc - mean(acc); % 3. 结构参数 Tn 1.0; % 自振周期秒矩形框架结构常见取值范围 zeta 0.05; % 阻尼比 % 4. 时域杜哈梅积分 [u_time, ~, ~] duhamel_response(t, acc_detrend, Tn, zeta); % 5. 频域求解 u_freq fft_response(t, acc_detrend, Tn, zeta); % 6. 对比图 figure(Color, w, Position, [100 100 800 500]); plot(t, u_time * 1000, b-, LineWidth, 1.2, DisplayName, Duhamel); hold on; plot(t, u_freq * 1000, r--, LineWidth, 1.2, DisplayName, FFT); xlabel(Time (s)); ylabel(Displacement (mm)); legend(Location, best); grid on; title(sprintf(Tn%.2fs, zeta%.2f, PGA%.2f g, Tn, zeta, max(abs(acc_detrend))/9.80665)); % 7. 峰值误差计算 peak_time max(abs(u_time)); peak_freq max(abs(u_freq)); err abs(peak_time - peak_freq) / peak_time * 100; fprintf(峰值位移: Duhamel%.4f mm, FFT%.4f mm, 误差%.2f%%\n, ... peak_time*1000, peak_freq*1000, err);主脚本里位数乘以1000是为了把米显示成毫米便于观察曲线幅值。PGA的计算用max(abs(acc_detrend))除以9.80665换算回g这是报告中常见的表述方式。两个结果误差小于2%就不用再排查了大于这个数就要回头调试。每次修改参数后重跑主脚本对比图会实时更新这是调试效率最高的方式。5.2 多周期扫描一次跑出T0.2秒到3秒的响应谱单一周期下的位移响应只是起步实际工程中更常见的是算一条位移反应谱也就是把结构自振周期从0.2秒扫到3秒记录每个周期的最大位移响应画成周期-位移曲线。这个过程同样基于前面两个函数只是外层包循环。位移反应谱对结构选型和抗震性能评估都有直接参考价值也是毕业设计里常出现的插图。%% 位移反应谱计算 T_list 0.2:0.1:3.0; spec_time zeros(size(T_list)); spec_freq zeros(size(T_list)); for i 1:length(T_list) [u_t, ~, ~] duhamel_response(t, acc_detrend, T_list(i), zeta); spec_time(i) max(abs(u_t)); u_f fft_response(t, acc_detrend, T_list(i), zeta); spec_freq(i) max(abs(u_f)); end figure(Color, w); plot(T_list, spec_time*1000, b-o, DisplayName, Duhamel); hold on; plot(T_list, spec_freq*1000, r-s, DisplayName, FFT); xlabel(Natural period Tn (s)); ylabel(Spectral displacement (mm)); legend(Location, northwest); grid on;这段循环在Matlab里跑起来很快3000点数据乘30个周期也就是几秒钟的事。如果你发现spec_time和spec_freq两条线在某个周期区间出现了系统性偏差比如短周期段FFT结果一直偏大检查一下是不是FFT的频率分辨率不够。频率分辨率等于Fs/n当n1500、Fs100时分辨率是0.067Hz相比结构自振频率来说偏粗可能让频响函数在峰值附近的采样不够精确导致FFT结果偏小。这时把FFT点数从n补零到n*4分辨率提升4倍就能改善。5.3 为什么不建议直接调用matlab内置函数解决一切Matlab内置了lsim、lsim相关控制工具箱函数对线性系统求解也很快但我的建议是核心计算用自己的函数做。理由不是内置函数不好而是你在做科研或课程作业时需要输出算法公式对应的代码逻辑而内置函数把物理参数隐藏在状态空间矩阵里一旦公式推导出错很难定位。自己写递推解的过程中你被迫把结构动力学里的每个变量含义梳理一遍这对理解duhamel积分本身有好处。另外在毕业设计和论文写作中评审人更愿意看到你自己实现的算法代码加验证图而不是只调Simulink或者lsim的输出截图。自写代码还能方便地扩展到非线性体系到那时可以替换刚度项和阻尼项做逐步积分。这也是一套代码从课程作业走向科研工具的自然路径。6. 画图配色、rms对比与快速验证的收尾技巧6.1 用rms值而不是只看峰值对比峰值误差小于2%算通过但峰值只代表一个点上的信息两条曲线可能在中间某段时间有系统性偏差却不体现在峰值上。更严格的验证方式是计算整条时程的均方根误差rms(u_time - u_freq) / rms(u_time)。这个比值反映整体波形的一致性一般小于5%就说明时域和频域结果在工程精度水平上一致。把这个值打印在主脚本里比单看峰值更可靠。rms_err rms(u_time - u_freq) / rms(u_time) * 100; fprintf(RMS相对误差: %.2f%%\n, rms_err);如果rms误差大而峰值误差小先怀疑是不是记录末尾段两者在小幅值区出现相位差。相位差来自阻尼项公式的符号书写不一致检查两个函数里阻尼项的写法是否为2i*zeta*wn*w。另一种可能是fft_response里没有做fftshift处理导致频域乘法时频率轴和序列不匹配。Matlab不接受你直接给负频率段乘传递函数时位置错位所以务必保持频率轴从0开始线性递增的方式这在此前的代码中已经做对了。6.2 一个只用10秒的快速验证方法每次写完代码改完参数后我建议先跑一个模拟信号快速验证生成一个正弦扫频信号代替真实地震波频率从0.5Hz扫描到10Hz持续时长20秒然后丢进两个函数里做计算。理论上这个信号的解析位移响应没有闭合解但你可以做一个荒唐但有效的验证把结构自振周期设成0.1秒让结构近乎刚体此时位移响应应该近似等于地面位移的负值。地面位移可以通过对加速度积分两次得到虽然带漂移但起始段的局部形态应该能对上。如果连这种极限测试都过不了跑真实地震记录也一定是错的。t_sim (0:1999) * 0.01; f0 0.5; f1 10; acc_sim sin(2*pi*(f0*t_sim (f1-f0)/2 * t_sim.^2)); [T_sim, ~, ~] duhamel_response(t_sim, acc_sim, 0.1, 0.02);这段代码十秒钟内就能算完你可以用它做断点调试。如果T_sim曲线形态和预期的冲击响应接近再换真实地震记录。调试顺序是模拟信号先行真实数据后行避免把真实数据的预处理噪声和算法bug混在一起排查。6.3 图片输出与报告整理的规范最终出图时把figure的Color设为white用exportgraphics导出成PDF或PNG分辨率设置成300dpi以上确保论文插图清晰。图注要写清楚记录名称、震级、结构参数和PGA值不要只写“位移时程对比”这种空泛标题。另一张频谱图建议用loglog坐标系因为地震波加速度谱的纵轴动态范围大线性坐标下小峰值会被压平。这两张图组合起来基本就是一段完整的“数据-处理-验证”的证据链放进答辩PPT或实验报告里都有说服力。最后提醒一句AT2文件里的记录元数据要单独存好PEER数据库的RSN编号、震级、距离、Vs30这些参数在写报告时全部要用到别等画完图再回头翻网页。把读取函数返回值扩展成结构体里面塞上震级、台站和Vs30字段后续做数据筛选会省很多事。本文还有配套的精品资源点击获取