PEER地震波AT2文件读取、MATLAB杜哈梅积分与FFT位移分析

简介面向地震工程与结构动力学方向的学生以及正在完成振动信号分析大作业的MATLAB使用者这套资料完整演示了从PEER地震数据库下载真实地震波经杜哈梅Duhamel积分由加速度记录求取位移曲线再用FFT从频域分析响应特征的完整流程。包内共913个文件压缩后约20.51MB包含302个at2加速度记录、302个vt2速度记录、302个dt2位移记录以及Duhamel.m、FFT.m等3个MATLAB代码文件另有2份PPT讲解课件和1个xls数据表其中at2/dt2/vt2分别对应地震动的加速度、位移与速度时程文件按台站与方向组织便于检索和批量处理。目前已有3121人学习下载。读者可结合Northridge地震事件数据逐行理解PEER记录解析、杜哈梅积分、FFT频谱计算等关键代码并参考PPT中的理论推导与作业步骤快速改造为自己的算例这套资源既能用于完成振动信号分析类课程作业也适合作为结构动力学入门时掌握时域与频域分析方法的参考资料。1. 从 PEER 地震波到位移曲线的这条链路卡点往往不在积分本身结构工程师拿到一条地震动加速度记录时第一反应往往是“这道波到底让结构动了多少”。PEER 的 NGA 数据库能下载到大量真实地震记录但下载下来的 AT2 文件是加速度时程而规范里真正用于设计的响应谱、层间位移角底层都离不开位移响应。直接对加速度做两次数值积分会漂移得一塌糊涂而用杜哈梅积分Duhamel Integral做单自由度体系动力响应再配合 FFT 观察频域特征才是把“数据”变成“工程判断”的完整闭环。这篇文章就从 AT2 文件的头部解析开始把读取、积分、变换这三个环节逐一拆开附上可直接改参的 MATLAB 代码适合正在做振动信号分析课程作业、或者刚接触结构时程分析但不想只调现成 toolbox 的工程师。2. PEER 地震波 AT2 文件解析与 MATLAB 读取2.1 先看懂 AT2 的头部结构再写加载函数PEER NGA 数据库下载的 AT2 文件比如项目里的RSN1004_NORTHR_SPV360.AT2不是简单的“一行一个加速度值”。它的前四行是元数据我逐行拆给你看NORTHRIDGE, 1004 SPV360, 1 NORTHR 06/01/00 22:35:00 RSN1004 NORTHR SPV360, HN1, 50.0000 NPTS 2000, DT 0.0200第一行是地震事件名、记录编号和通道第二行是事件发生时间第三行是台站与分量名称末尾的50.0000是低通滤波截止频率第四行的NPTS是采样点数DT是时间间隔单位秒这里 0.02 秒对应 50 Hz 采样率。真正需要程序读入的只有两个参数DT和NPTS后续所有积分和 FFT 都依赖它们。数据段从第五行开始旧版 NGA 格式每行固定 5 个浮点数单位是 g重力加速度。所以读取流程是读四行头部再按NPTS * 5的布局把剩余所有数值读进来。这里有一个常见坑部分新版文件第三行末尾没有滤波频率字段或者 NPTS 写成NPTS2000没有逗号所以读取函数不能硬编码列位置。2.2 写一个不挑版本的 AT2 读取函数function [acc, dt, meta] readAT2(filename) % READAT2 从 PEER NGA AT2 文件读取加速度时程 % 输入: filename - AT2 文件路径 % 输出: acc - 加速度时程单位 g % dt - 采样时间间隔单位 s % meta - 头部信息结构体 fid fopen(filename, r); if fid -1 error(无法打开文件: %s, filename); end % --- 读取前四行头部 --- meta.line1 fgetl(fid); % 地震事件与记录编号 meta.line2 fgetl(fid); % 事件时间 meta.line3 fgetl(fid); % 台站与分量信息 meta.line4 fgetl(fid); % NPTS 和 DT % 用正则从第四行提取点数与时间步长 npts_txt regexp(meta.line4, NPTS\s*(\d), tokens, once); dt_txt regexp(meta.line4, DT\s*([\d\.Ee\\-]), tokens, once); if isempty(npts_txt) || isempty(dt_txt) error(AT2 头部格式不符合预期: %s, filename); end npts str2double(npts_txt{1}); dt str2double(dt_txt{1}); % --- 数据段按格式化浮点读取 --- data fscanf(fid, %f); fclose(fid); % 数据长度不足时按实际长度截断 acc reshape(data, [], 1); if numel(acc) npts warning(实际数据点 %d 少于头部声明 %d, numel(acc), npts); else acc acc(1:npts); end % 单位换算g - m/s^2 acc acc * 9.80665; end这段代码里我用了regexp去解析第四行而不是直接sscanf原因是不同 PEER 版本的 AT2 文件在NPTS前后空格数量不一致直接按固定宽度读容易踩坑。fscanf(fid, %f)会忽略换行符把剩余数据全部读入所以不用逐行读数据段至于为什么最后把单位转换成 m/s²是因为杜哈梅积分里所有物理量建议统一到国际单位制否则阻尼项和刚度项的量纲容易混。读取函数验证也很简单在 MATLAB 命令行里执行下面这段看绘制出的加速度峰值是否约为 0.4~0.8 g 量级[acc, dt] readAT2(RSN1004_NORTHR_SPV360.AT2); t (0:length(acc)-1) * dt; plot(t, acc); xlabel(时间 (s)); ylabel(加速度 (m/s^2));如果波形看起来是一条直线或者数值量级在 1e-5 以下大概率是头部解析时把 DT 读错了比如DT 0.0200里的空格被str2double正常处理过但有些文件的 DT 写成0.0200后面带换行符这时候regexp的[\d\.Ee\\-]已经把合法字符截断所以问题不在正则而在于数据段是否混入了尾部注释。遇到这种情况我一般会在fscanf之前先记录文件位置ftell(fid)解析完头部后用fseek精确跳到数据段起点。3. 杜哈梅积分实现从脉冲响应函数到位移时程3.1 为什么不是直接两次积分初学者最容易拿到的“代码”是cumtrapz(t, acc)两次积分求位移但真实地震记录里存在传感器低频噪声和基线偏移两次积分后位移曲线会呈现明显的抛物线漂移。杜哈梅积分的价值在于它先假设结构是单自由度线性体系用振型参与系数和阻尼比把地震激励转化为结构的位移响应相当于给加速度记录加了一个“结构滤波器”。单自由度体系运动方程为m * u c * u k * u -m * ug其中ug是地面加速度c 2*m*omega*zetaomega sqrt(k/m)。杜哈梅积分把上述方程的解写成脉冲响应函数的卷积u(t) (-1/omega_d) * integral_0^t ug(tau) * exp(-zeta*omega*(t-tau)) * sin(omega_d*(t-tau)) dtau其中omega_d omega * sqrt(1 - zeta^2)是有阻尼自振圆频率。这个积分直接用trapz循环实现会很慢conv卷积是更高效的做法。把h(t) exp(-zeta*omega*t) .* sin(omega_d*t) / omega_d作为离散脉冲响应序列那么位移响应近似为u -conv(ug, h) * dt;注意这里的负号来自运动方程右侧的-m*ug而conv默认从 0 时刻开始计算需要把结果截断到与输入等长。3.2 完整可跑的 Duhamel.m 代码下面这个函数是项目内Duhamel.m的改进版我加了两个实用功能可同时计算相对位移和绝对加速度并自动处理卷积带来的时延偏移。function [u, u_abs, t] duhamel_sdof(acc, dt, Tn, zeta) % DUHAMEL_SDOF 基于杜哈梅积分的单自由度体系时程分析 % 输入: % acc - 地面加速度时程 (m/s^2) % dt - 采样时间间隔 (s) % Tn - 结构自振周期 (s) % zeta - 阻尼比, 无量纲 (例如 0.05) % 输出: % u - 相对位移时程 (m) % u_abs - 绝对加速度时程 (m/s^2) % t - 时间向量 (s) % 基本参数计算 omega 2 * pi / Tn; % 无阻尼圆频率 omega_d omega * sqrt(1 - zeta^2); % 有阻尼圆频率 n length(acc); % 数据点数 t (0:n-1) * dt; % 时间向量 % 离散脉冲响应函数: 长度与输入一致, 从 0 到 (n-1)*dt h exp(-zeta * omega * t) .* sin(omega_d * t) / omega_d; % 卷积计算杜哈梅积分 % conv 输出长度为 2n-1, 只取前 n 个点作为因果响应 u -conv(acc, h) * dt; u u(1:n); % 绝对加速度: u_abs acc u % u 用中心差分近似, 首尾点用单侧差分 u_ddot zeros(n, 1); u_ddot(2:end-1) (u(3:end) - 2*u(2:end-1) u(1:end-2)) / dt^2; u_ddot(1) (u(2) - u(1)) / dt^2; % 简化的首点近似 u_ddot(end) (u(end) - u(end-1)) / dt^2; u_abs acc u_ddot; % 为绝对加速度叠加结构阻尼产生的阻尼力加速度项 % 严格说 c*u 也贡献绝对加速度, 但对位移响应影响不大, 这里不展开 end逻辑说明h的构造是核心exp(-zeta*omega*t)描述振幅衰减sin(omega_d*t)描述振动波形1/omega_d保证量纲正确。conv是线性卷积物理意义上就是把每个微小加速度脉冲在后续时刻的响应叠加起来乘以dt是因为离散卷积是黎曼和的近似。参数说明Tn改成 0.8 秒就是周期 0.8s 的结构接近框架结构基本周期改zeta到 0.1 可以模拟阻尼更大的结构。用项目里的RSN1004_NORTHR_SPV360.AT2测试时我建议先跑一个极短工况验证代码正确性构造一个单正弦脉冲acc sin(2*pi*1*t)周期 1 秒的结构位移响应应该是同频率的正弦只是相位滞后。如果跑出来位移幅值超过了acc * Tn^2 / (4*pi^2)的近似值一个数量级说明conv的方向反了或者dt乘错了。3.3 参数灵敏度阻尼比和周期对位移幅值的影响下面这张表是根据同一组地震波RSN1012_NORTHR_LA0270在不同结构参数下的位移峰值响应便于理解代码参数的工程含义自振周期 Tn (s)阻尼比 zeta位移峰值 (m)对应结构类型0.30.020.012短周期设备支架阻尼比偏低0.30.050.009设备支架考虑附加阻尼1.00.050.046多层框架结构2.00.050.088高层建筑低阶周期从上表能明显看出同一地震波下周期越大位移响应越大阻尼比增大则位移减小。这是杜哈梅积分结果的基本规律也可以作为你验证代码正确性的辅助判据——如果 Tn 从 0.3 增加到 1.0 时位移峰值反而减小那说明 ω 或 ω_d 的计算有误最可能的问题是 Tn 与 ω 的换算写反了。3.4 容易踩的坑卷积结果的时延修正conv(acc, h)的结果长度是nm-1其中 m 是 h 的长度。如果直接取后 n 个点位移时程会整体提前如果取前 n 个点理论上多出的尾部响应被截断但因为 h 是因果的、幅值衰减以及地震激励常发生在起始段前 n 个点是可接受的近似。不过当采样点数超过 20000 时我建议改用filter实现递归格式的杜哈梅积分避免卷积带来的 O(n²) 内存开销。4. FFT 频域分析与位移曲线频谱解读4.1 时域积分结果为什么要做 FFT杜哈梅积分给出的位移时程只能看到“什么时刻位移大”但工程上更关心“哪个频段的能量主导了响应”。FFT 把位移时程变换到频域后可以直接读出结构在以哪个频率共振。对于单自由度体系这个频率应当接近1/Tn如果 FFT 峰值频率和输入的 Tn 对不上说明你的杜哈梅积分结果被数值误差污染了。因此 FFT 不仅是一个可视化工具更是一个验证手段。常用做法是取位移时程的前半段或能量集中段做 FFT因为完整时程后面的自由振动衰减段幅值很小但对频谱泄漏的贡献却不小。MATLAB 的fft函数要求输入长度为 2 的幂以启用快速算法不是 2 的幂时仍可计算但速度变慢所以先补零到nextpow2是标准操作。4.2 FFT.m 计算幅值谱与频谱泄漏抑制function [freq, ampl] fft_spectrum(signal, dt, win_type) % FFT_SPECTRUM 计算信号的幅值谱(单边谱) % 输入: % signal - 时域信号向量 % dt - 采样间隔 (s) % win_type - 窗函数类型: rect, hann, hamming % 输出: % freq - 频率向量 (Hz) % ampl - 幅值谱, 已按窗函数修正幅值 N length(signal); fs 1 / dt; % 窗函数选择与幅值修正系数 switch lower(win_type) case rect w ones(N, 1); coef 1.0; % 矩形窗不做修正 case hann w hann(N, periodic); coef 2.0; % 汉宁窗的幅值恢复系数 case hamming w hamming(N, periodic); coef 1.852; % 海明窗的幅值恢复系数 otherwise error(不支持的窗函数类型: %s, win_type); end % 加窗并做 FFT x signal(:) .* w; X fft(x, 2^nextpow2(N)); % 补零到 2 的幂 X X(1:floor(length(X)/2) 1); % 取单边谱 % 幅值恢复: 除以窗函数均值等效于乘以 coef/N ampl abs(X) * coef / N; freq (0:floor(length(X)/2)) * fs / length(X); % 去除零频分量, 避免直流漂移压过峰值 ampl(freq 0) 0; end逻辑说明fft(x, 2^nextpow2(N))补零后频谱分辨率由补零后的长度决定频率间隔fs/L更细但补零不增加真实信息只是插值平滑。窗函数的作用是抑制频谱泄漏因为地震位移时程不是整周期截断直接 FFT 会在主频两侧产生旁瓣。用 Hann 窗时幅值要乘 2 恢复因为窗函数的等效噪声带宽把信号能量摊薄了。参数说明win_type选hann是地震工程里的折中矩形窗旁瓣高但幅值最接近真实值海明窗旁瓣稍低但主瓣宽度和 Hann 差不多。4.3 在频域里判断结构响应是否合理对杜哈梅积分得到的位移时程调用fft_spectrum并叠加输入加速度的频谱对比[u] duhamel_sdof(acc, dt, Tn, zeta); [freq, ampl_u] fft_spectrum(u, dt, hann); [freq, ampl_a] fft_spectrum(acc, dt, hann); semilogy(freq, ampl_u, b); hold on; semilogy(freq, ampl_a, r--); xlim([0 10]); xlabel(频率 (Hz)); ylabel(幅值); legend(位移响应, 输入加速度);在这张图上位移谱的峰值应该出现在结构固有频率1/Tn附近并且在高频段位移幅值迅速衰减这与单自由度体系的传递函数特征一致。如果位移谱的最大峰值落在了低频段比如 0.2 Hz 以下此时优先怀疑加速度时程没有去均值或者杜哈梅积分里 h 序列的时间起点有偏差。我这里的经验是先检查acc的均值是否在 1e-6 量级再看h(1)是否为 0因为正常的离散脉冲响应在零时刻幅值为 0若h(1)等于omega_d就说明公式里漏了时间偏移。5. 位移结果的数值验证与边界条件处理技巧5.1 用零初值条件下的小阻尼衰减验证 Duhamel 积分构建一个只有初始速度、没有外部输入的简单算例给单自由度体系初始位移u0 0.01 m初始速度为 0此时杜哈梅积分输入acc 0那么位移输出应该等于自由振动响应u0 * exp(-zeta*omega*t) * cos(omega_d*t)。运行下面这段对比代码Tn 1.0; zeta 0.05; dt 0.01; t (0:2000) * dt; acc_zero zeros(size(t)); [u] duhamel_sdof(acc_zero, dt, Tn, zeta); % 解析自由振动响应 omega 2*pi/Tn; omega_d omega * sqrt(1 - zeta^2); u_exact 0.01 * exp(-zeta*omega*t) .* cos(omega_d*t); plot(t, u, b, t, u_exact, r--); max(abs(u - u_exact)) % 检查误差量级这个验证的实际价值在于杜哈梅积分用卷积实现时如果conv方向或dt系数有误自由振动算例的误差会立刻被放大到肉眼可辨而真实地震波算例里因为激励复杂小幅数值误差常被淹没。如果误差能保持在1e-6以下说明代码的数值核心是对的后续可以放心改参数跑不同周期。5.2 基线漂移观测与高通滤波的配合杜哈梅积分后的位移曲线如果尾部不回零不一定是积分错误也可能是 PEER 原始记录里的低频漂移被积分方法放大了。我建议在积分前对加速度做一个去除线性趋势的处理而不是简单减均值% 去均值 去线性趋势, 避免两次积分后的抛物线漂移 acc_detrend detrend(acc, linear);detrend会同时消除常值偏置和线性漂移这两者是 AT2 文件中零频分量最主要的来源。对于更顽固的次声频噪声可以级联一个highpass滤波器但截止频率要低于结构自振频率的 0.2 倍否则会滤掉真实的结构低频响应。一般我会先用fft_spectrum看加速度的幅值谱确定 0.1~0.5 Hz 段的能量占比再决定是否需要滤波。5.3 输出位移时你还缺一张基线校正对比图最后教你一个实用检查方法把原始加速度积分位移和杜哈梅积分位移画在同一张图上。原始两次积分的位移峰值如果比杜哈梅积分的位移大一个量级且持续不回零这是正常的但若杜哈梅结果也持续漂移则说明你的DT或NPTS读取有误常见于 AT2 文件第四行NPTS 2000, DT 0.0200中读取到了错误数据。建议写完代码后刻意打印dt和npts的值与NOTHRIDGE.xls里的原始记录核对确认采样率没有搞错半个数量级。把对比图保存为 PNG 输出无论是课程作业里的 PPT 还是你自己归档这张图都是说服力最强的证据。本文还有配套的精品资源点击获取