简介基于MATLAB的ERA环境激励模态识别资源包面向结构动力学与健康监测领域的研究者和工程师解决从随机环境响应中提取固有频率、阻尼比和振型等模态参数的工程问题。资源以ERA工具箱为核心包含ERA主程序、多类激励数据生成脚本、模态质量与MAC计算模块等6个m功能文件配套3个asv自动保存版本及12个txt说明文档共21个文件压缩后仅18KB轻量易用。目前已有700余人学习下载。使用者可依据说明快速运行演示程序结合普通、阶跃、正弦等数据生成器模拟不同环境激励再通过模态参数、模态质量及MAC等输出结果验证识别效果适合需要快速开展模态参数识别实验或教学演示的群体。1. ERA 模态识别环境激励下真正能落地的 MATLAB 工具箱拿到 ERA.rar 这套源码时先纠正一个常见的认知偏差ERA 并不是摘要里写的 Empirical Mode Decomposition那叫 EMD而是 Eigensystem Realization Algorithm特征系统实现算法。这是 NASA 在 1980 年代为大型空间结构模态辨识开发的时域方法核心思想是从自由衰减响应或脉冲响应中构造 Hankel 矩阵用 SVD 做最小实现再从状态矩阵的特征分解里直接读出频率、阻尼和振型。相比频域峰值拾取法PP和频域分解法FDDERA 不依赖频响函数估计对环境激励下的弱激励、密集模态、高阻尼结构都更稳。这套 MATLAB 工具箱把从激励数据生成、ERA 求解到 MAC 验证的完整链路都打包了适合正在做结构健康监测、桥梁/风机/建筑现场实测的工程师也适合研究生拿真实信号练手。2. ERA 的数学骨架从 Hankel 矩阵到最小实现2.1 为什么 ERA 适合环境激励数据环境激励风、交通、地脉动的特点是输入不可测、幅值随机、频带较宽。传统方法需要已知输入做频响函数FRF估计环境激励下只能靠输出数据。ERA 走的是另一条路先通过自然激励技术NExT把响应的互相关函数当作自由衰减响应再喂给 ERA 做状态空间辨识。工具箱里的 norsindata_generater.m 生成的就是这种随机响应数据配合 normaldata_generater.m 的高斯白噪声激励正好模拟了环境激励的场景。2.2 核心推导Hankel 矩阵与 SVD 截断ERA 的起点是脉冲响应序列 (h(k))。对多输入多输出系统构造分块 Hankel 矩阵[ H(k-1) \begin{bmatrix} h(k) h(k1) \cdots h(ks-1) \ h(k1) h(k2) \cdots h(ks) \ \vdots \vdots \ddots \vdots \ h(kr-1) h(kr) \cdots h(krs-2) \end{bmatrix} ]其中 (r) 是行块数(s) 是列块数都必须大于系统阶数的两倍。对 (H(0)) 做 SVD[ H(0) U \Sigma V^T ]用前 (n) 个奇异值截断(n) 为系统真实阶数的两倍得到最小实现的状态矩阵[ A \Sigma_n^{-1/2} U_n^T H(1) V_n \Sigma_n^{-1/2} ]对 (A) 做特征值分解特征值 (\lambda_i) 对应连续域极点[ s_i \frac{\ln(\lambda_i)}{\Delta t}, \quad f_i \frac{|s_i|}{2\pi}, \quad \zeta_i -\frac{\operatorname{Re}(s_i)}{|s_i|} ]2.3 在 MATLAB 里的等价实现工具箱的 ERA.m 核心逻辑可以简化为function [fn, zeta, phi] ERA(Y, fs, r, s, n) % Y: 脉冲响应或自由衰减响应维度为 no×nt % fs: 采样频率r/s: Hankel矩阵行/列块数n: 截断阶数 [no, nt] size(Y); H0 zeros(no*r, no*s); H1 zeros(no*r, no*s); for k 1:r for j 1:s H0((k-1)*no1:k*no, (j-1)*no1:j*no) Y(:, kj-1); H1((k-1)*no1:k*no, (j-1)*no1:j*no) Y(:, kj); end end [U, S, V] svd(H0, econ); Un U(:, 1:n); Sn S(1:n, 1:n); Vn V(:, 1:n); A sqrt(Sn) \ (Un * H1 * Vn) / sqrt(Sn); [Psi, D] eig(A); lambda diag(D); dt 1/fs; s_cont log(lambda) / dt; fn abs(s_cont) / (2*pi); zeta -real(s_cont) ./ abs(s_cont); C Y(:, 1:r) * Vn * sqrt(Sn) * Psi; % 输出矩阵C phi C; % 振型取输出矩阵列向量 end代码逻辑说明双重循环构造 H0 和 H1 两个 Hankel 矩阵H1 相对 H0 整体后移一个时间步这是实现最小实现求解的关键。svd(H0, econ)做经济型 SVD省内存当no*r远大于 200 阶时优势很明显。sqrt(Sn) \ (Un * H1 * Vn) / sqrt(Sn)等价于 (\Sigma_n^{-1/2} U_n^T H(1) V_n \Sigma_n^{-1/2})用左除避免显式求逆。特征值取对数转连续域时注意虚部可能跨 (\pm\pi) 边界工程上一般只取正实部对应的模态。3. 工具箱结构与调用链拆解3.1 文件清单与职责解压 ERA.rar 后核心文件和职责如下表文件职责ERA.m主算法函数输入响应/脉冲序列输出频率、阻尼、振型ERA_StartDemo.m演示入口串起数据生成-识别-后处理全过程normaldata_generater.m生成高斯白噪声激励下的响应数据sindata_generater.m生成正弦激励下的稳态响应数据norsindata_generater.m生成随机正弦混合激励响应stepdata_generater.m生成阶跃激励响应MAC.txt / MAC.m模态置信准则计算验证振型一致性System_Modal_Parameters.txt理论模态参数真值用于对比System_Eigenvalue.txt理论特征值System_Mass_M.txt / Modal_Mass_M.txt质量矩阵和模态质量sampling_frequencyHz.txt采样频率YData.txt实测/仿真响应数据ERA_readme.txt使用说明注意一下System_Modal_Parameters.txt 存的是结构理论值不是 ERA 识别结果。调试时拿它做真值对比识别偏差。3.2 ERA_StartDemo.m 的调用流程% ERA_StartDemo.m 关键段落 % 加载仿真数据 Y load(YData.txt); % 响应矩阵每行一个测点 fs load(sampling_frequencyHz.txt); % 生成脉冲响应ERA 需要自由衰减或脉冲响应 h norsindata_generater(Y, fs); % 混合激励下的响应预处理 % 设置 Hankel 矩阵维度 r 20; % 行块数 s 20; % 列块数 n 6; % 截断阶数2倍模态数 % 执行 ERA [fn, zeta, phi] ERA(h, fs, r, s, n); % MAC 验证 mac MAC(phi, phi_ref); % phi_ref 从 System_Shape.txt 读取逻辑说明norsindata_generater在这里不是简单的数据生成器它内部做了互相关处理输出近似脉冲响应序列这是环境激励下应用 ERA 的前提。r和s取 20 是经验值对 3 阶主模态的系统阶数 6 的 3 倍以上就够。取值太小2n会让 Hankel 矩阵病态太大50会引入噪声子空间。n取实际模态数的 2 倍因为状态空间方程每个模态对应一对共轭复极点。4. 跑通识别流程与参数调优4.1 用 demo 脚本验证全链路在 MATLAB R2018b 及以上版本直接运行 ERA_StartDemo.m观察输出。以三自由度弹簧-质量系统为例工具箱内置算例采样频率设为 200 Hz理论频率通常在 5-20 Hz 区间。运行后应看到识别模态1: 5.02 Hz, 阻尼比 0.80% 理论模态1: 5.00 Hz, 阻尼比 0.80% 识别模态2: 12.11 Hz, 阻尼比 1.20% 理论模态2: 12.07 Hz, 阻尼比 1.20% 识别模态3: 18.85 Hz, 阻尼比 1.50%如果识别频率和理论频率偏差超过 2%优先检查sampling_frequencyHz.txt是否和生成数据时一致。工具箱的 normaldata_generater.m 在生成数据时会把采样频率写入该文件如果手动拼接数据源这是最容易出错的点。4.2 关键参数怎么定参数推荐范围判断依据采样频率 fs信号最高关注频率的 5~10 倍低于 5 倍阻尼比识别偏差大Hankel 行块 r2n5 ~ 3n过小丢失模态过大增加噪声Hankel 列块 s与 r 相同或略大保证 H0 接近方阵SVD 稳定截断阶数 n2×(关注模态数)若奇异值谱无显著跌落逐步加 2 试探数据长度至少包含 20 个最低频周期环境激励下建议 50 个周期以上调试建议先扫 n 从 2 到 20 的步进看频率估计是否收敛。真正可信的模态在 n 变化时频率漂移应小于 0.5%。这是区分真实模态和噪声模态的常用手段。5. 从 NExT 到稳态图环境激励实操的关键步骤5.1 NExT 与 ERA 的联系环境激励下没有脉冲响应可用必须先从实测响应计算互相关。NExTNatural Excitation Technique的结论是在白噪声激励下两个测点响应的互相关函数满足齐次运动方程可以当作自由衰减响应来用。工具箱的 metode就是先做互相关再进 ERA。实操中我一般用参考点法function [R, tau] NExT_crosscorr(Y, ref_idx, max_lag) % Y: 响应矩阵 [测点 × 时间] % ref_idx: 参考测点序号选振幅大、信噪比高的点 [n_pts, n_t] size(Y); R zeros(n_pts, max_lag); for i 1:n_pts [R(i,:), tau] xcorr(Y(i,:), Y(ref_idx,:), max_lag, biased); end R R(:, max_lag1:end); % 取正滞后部分 end参考点选不好会直接丢模态。经验规则参考点要落在所有关注模态振型的非节点位置否则该模态在互相关里能量太低ERA 识别不出来。多试几个候选点看哪组 MAC 值整体最高。5.2 稳态图的实现技巧单次 ERA 识别结果未必可信工程上常用稳态图筛选不断增加截断阶数 n把识别出的频率画在同一张图上真实模态会在某个 n 值后趋于稳定噪声模态则四处漂移。以下这段适用于批量扫参n_list 2:2:30; f_stable []; for n n_list [fn, zeta, phi] ERA(h, fs, 20, 20, n); f_stable [f_stable; fn(:)]; % 稳定性判据频率变化 1%阻尼变化 5% end绘图后能看见 3-4 条竖直的稳定线这些就是物理模态斜线是噪声极点直接忽略。增加了阻尼比的稳定性判据后噪点大幅减少但要注意阻尼比估计本身方差大判据放宽到 5% 比较实用。5.3 MAC 校验和阻尼修正的坑MAC.txt 存的是模态置信准则矩阵计算公式是[ MAC(i,j) \frac{|\phi_i^T \phi_j|^2}{(\phi_i^T \phi_i)(\phi_j^T \phi_j)} ]对角元接近 1、非对角元接近 0 说明振型分离得干净。如果两个模态的 MAC 值超过 0.7基本上可以判定识别出了重复模态或振型混叠这时需要增加测点或重新选参考点。环境激励识别的阻尼比通常偏大原因是互相关函数尾部信噪比低ERA 拟合时会把噪声衰减也算进去。对策是把互相关函数截断到衰减至峰值 10% 左右的长度再喂给 ERA能修正约 0.3~0.5 个百分点的阻尼偏差。MATLAB 里直接取 R 序列前 1/3 段即可代价是频率分辨率略有下降。本文还有配套的精品资源点击获取 ERA模态识别实战：环境激励下基于MATLAB的特征系统实现算法
简介基于MATLAB的ERA环境激励模态识别资源包面向结构动力学与健康监测领域的研究者和工程师解决从随机环境响应中提取固有频率、阻尼比和振型等模态参数的工程问题。资源以ERA工具箱为核心包含ERA主程序、多类激励数据生成脚本、模态质量与MAC计算模块等6个m功能文件配套3个asv自动保存版本及12个txt说明文档共21个文件压缩后仅18KB轻量易用。目前已有700余人学习下载。使用者可依据说明快速运行演示程序结合普通、阶跃、正弦等数据生成器模拟不同环境激励再通过模态参数、模态质量及MAC等输出结果验证识别效果适合需要快速开展模态参数识别实验或教学演示的群体。1. ERA 模态识别环境激励下真正能落地的 MATLAB 工具箱拿到 ERA.rar 这套源码时先纠正一个常见的认知偏差ERA 并不是摘要里写的 Empirical Mode Decomposition那叫 EMD而是 Eigensystem Realization Algorithm特征系统实现算法。这是 NASA 在 1980 年代为大型空间结构模态辨识开发的时域方法核心思想是从自由衰减响应或脉冲响应中构造 Hankel 矩阵用 SVD 做最小实现再从状态矩阵的特征分解里直接读出频率、阻尼和振型。相比频域峰值拾取法PP和频域分解法FDDERA 不依赖频响函数估计对环境激励下的弱激励、密集模态、高阻尼结构都更稳。这套 MATLAB 工具箱把从激励数据生成、ERA 求解到 MAC 验证的完整链路都打包了适合正在做结构健康监测、桥梁/风机/建筑现场实测的工程师也适合研究生拿真实信号练手。2. ERA 的数学骨架从 Hankel 矩阵到最小实现2.1 为什么 ERA 适合环境激励数据环境激励风、交通、地脉动的特点是输入不可测、幅值随机、频带较宽。传统方法需要已知输入做频响函数FRF估计环境激励下只能靠输出数据。ERA 走的是另一条路先通过自然激励技术NExT把响应的互相关函数当作自由衰减响应再喂给 ERA 做状态空间辨识。工具箱里的 norsindata_generater.m 生成的就是这种随机响应数据配合 normaldata_generater.m 的高斯白噪声激励正好模拟了环境激励的场景。2.2 核心推导Hankel 矩阵与 SVD 截断ERA 的起点是脉冲响应序列 (h(k))。对多输入多输出系统构造分块 Hankel 矩阵[ H(k-1) \begin{bmatrix} h(k) h(k1) \cdots h(ks-1) \ h(k1) h(k2) \cdots h(ks) \ \vdots \vdots \ddots \vdots \ h(kr-1) h(kr) \cdots h(krs-2) \end{bmatrix} ]其中 (r) 是行块数(s) 是列块数都必须大于系统阶数的两倍。对 (H(0)) 做 SVD[ H(0) U \Sigma V^T ]用前 (n) 个奇异值截断(n) 为系统真实阶数的两倍得到最小实现的状态矩阵[ A \Sigma_n^{-1/2} U_n^T H(1) V_n \Sigma_n^{-1/2} ]对 (A) 做特征值分解特征值 (\lambda_i) 对应连续域极点[ s_i \frac{\ln(\lambda_i)}{\Delta t}, \quad f_i \frac{|s_i|}{2\pi}, \quad \zeta_i -\frac{\operatorname{Re}(s_i)}{|s_i|} ]2.3 在 MATLAB 里的等价实现工具箱的 ERA.m 核心逻辑可以简化为function [fn, zeta, phi] ERA(Y, fs, r, s, n) % Y: 脉冲响应或自由衰减响应维度为 no×nt % fs: 采样频率r/s: Hankel矩阵行/列块数n: 截断阶数 [no, nt] size(Y); H0 zeros(no*r, no*s); H1 zeros(no*r, no*s); for k 1:r for j 1:s H0((k-1)*no1:k*no, (j-1)*no1:j*no) Y(:, kj-1); H1((k-1)*no1:k*no, (j-1)*no1:j*no) Y(:, kj); end end [U, S, V] svd(H0, econ); Un U(:, 1:n); Sn S(1:n, 1:n); Vn V(:, 1:n); A sqrt(Sn) \ (Un * H1 * Vn) / sqrt(Sn); [Psi, D] eig(A); lambda diag(D); dt 1/fs; s_cont log(lambda) / dt; fn abs(s_cont) / (2*pi); zeta -real(s_cont) ./ abs(s_cont); C Y(:, 1:r) * Vn * sqrt(Sn) * Psi; % 输出矩阵C phi C; % 振型取输出矩阵列向量 end代码逻辑说明双重循环构造 H0 和 H1 两个 Hankel 矩阵H1 相对 H0 整体后移一个时间步这是实现最小实现求解的关键。svd(H0, econ)做经济型 SVD省内存当no*r远大于 200 阶时优势很明显。sqrt(Sn) \ (Un * H1 * Vn) / sqrt(Sn)等价于 (\Sigma_n^{-1/2} U_n^T H(1) V_n \Sigma_n^{-1/2})用左除避免显式求逆。特征值取对数转连续域时注意虚部可能跨 (\pm\pi) 边界工程上一般只取正实部对应的模态。3. 工具箱结构与调用链拆解3.1 文件清单与职责解压 ERA.rar 后核心文件和职责如下表文件职责ERA.m主算法函数输入响应/脉冲序列输出频率、阻尼、振型ERA_StartDemo.m演示入口串起数据生成-识别-后处理全过程normaldata_generater.m生成高斯白噪声激励下的响应数据sindata_generater.m生成正弦激励下的稳态响应数据norsindata_generater.m生成随机正弦混合激励响应stepdata_generater.m生成阶跃激励响应MAC.txt / MAC.m模态置信准则计算验证振型一致性System_Modal_Parameters.txt理论模态参数真值用于对比System_Eigenvalue.txt理论特征值System_Mass_M.txt / Modal_Mass_M.txt质量矩阵和模态质量sampling_frequencyHz.txt采样频率YData.txt实测/仿真响应数据ERA_readme.txt使用说明注意一下System_Modal_Parameters.txt 存的是结构理论值不是 ERA 识别结果。调试时拿它做真值对比识别偏差。3.2 ERA_StartDemo.m 的调用流程% ERA_StartDemo.m 关键段落 % 加载仿真数据 Y load(YData.txt); % 响应矩阵每行一个测点 fs load(sampling_frequencyHz.txt); % 生成脉冲响应ERA 需要自由衰减或脉冲响应 h norsindata_generater(Y, fs); % 混合激励下的响应预处理 % 设置 Hankel 矩阵维度 r 20; % 行块数 s 20; % 列块数 n 6; % 截断阶数2倍模态数 % 执行 ERA [fn, zeta, phi] ERA(h, fs, r, s, n); % MAC 验证 mac MAC(phi, phi_ref); % phi_ref 从 System_Shape.txt 读取逻辑说明norsindata_generater在这里不是简单的数据生成器它内部做了互相关处理输出近似脉冲响应序列这是环境激励下应用 ERA 的前提。r和s取 20 是经验值对 3 阶主模态的系统阶数 6 的 3 倍以上就够。取值太小2n会让 Hankel 矩阵病态太大50会引入噪声子空间。n取实际模态数的 2 倍因为状态空间方程每个模态对应一对共轭复极点。4. 跑通识别流程与参数调优4.1 用 demo 脚本验证全链路在 MATLAB R2018b 及以上版本直接运行 ERA_StartDemo.m观察输出。以三自由度弹簧-质量系统为例工具箱内置算例采样频率设为 200 Hz理论频率通常在 5-20 Hz 区间。运行后应看到识别模态1: 5.02 Hz, 阻尼比 0.80% 理论模态1: 5.00 Hz, 阻尼比 0.80% 识别模态2: 12.11 Hz, 阻尼比 1.20% 理论模态2: 12.07 Hz, 阻尼比 1.20% 识别模态3: 18.85 Hz, 阻尼比 1.50%如果识别频率和理论频率偏差超过 2%优先检查sampling_frequencyHz.txt是否和生成数据时一致。工具箱的 normaldata_generater.m 在生成数据时会把采样频率写入该文件如果手动拼接数据源这是最容易出错的点。4.2 关键参数怎么定参数推荐范围判断依据采样频率 fs信号最高关注频率的 5~10 倍低于 5 倍阻尼比识别偏差大Hankel 行块 r2n5 ~ 3n过小丢失模态过大增加噪声Hankel 列块 s与 r 相同或略大保证 H0 接近方阵SVD 稳定截断阶数 n2×(关注模态数)若奇异值谱无显著跌落逐步加 2 试探数据长度至少包含 20 个最低频周期环境激励下建议 50 个周期以上调试建议先扫 n 从 2 到 20 的步进看频率估计是否收敛。真正可信的模态在 n 变化时频率漂移应小于 0.5%。这是区分真实模态和噪声模态的常用手段。5. 从 NExT 到稳态图环境激励实操的关键步骤5.1 NExT 与 ERA 的联系环境激励下没有脉冲响应可用必须先从实测响应计算互相关。NExTNatural Excitation Technique的结论是在白噪声激励下两个测点响应的互相关函数满足齐次运动方程可以当作自由衰减响应来用。工具箱的 metode就是先做互相关再进 ERA。实操中我一般用参考点法function [R, tau] NExT_crosscorr(Y, ref_idx, max_lag) % Y: 响应矩阵 [测点 × 时间] % ref_idx: 参考测点序号选振幅大、信噪比高的点 [n_pts, n_t] size(Y); R zeros(n_pts, max_lag); for i 1:n_pts [R(i,:), tau] xcorr(Y(i,:), Y(ref_idx,:), max_lag, biased); end R R(:, max_lag1:end); % 取正滞后部分 end参考点选不好会直接丢模态。经验规则参考点要落在所有关注模态振型的非节点位置否则该模态在互相关里能量太低ERA 识别不出来。多试几个候选点看哪组 MAC 值整体最高。5.2 稳态图的实现技巧单次 ERA 识别结果未必可信工程上常用稳态图筛选不断增加截断阶数 n把识别出的频率画在同一张图上真实模态会在某个 n 值后趋于稳定噪声模态则四处漂移。以下这段适用于批量扫参n_list 2:2:30; f_stable []; for n n_list [fn, zeta, phi] ERA(h, fs, 20, 20, n); f_stable [f_stable; fn(:)]; % 稳定性判据频率变化 1%阻尼变化 5% end绘图后能看见 3-4 条竖直的稳定线这些就是物理模态斜线是噪声极点直接忽略。增加了阻尼比的稳定性判据后噪点大幅减少但要注意阻尼比估计本身方差大判据放宽到 5% 比较实用。5.3 MAC 校验和阻尼修正的坑MAC.txt 存的是模态置信准则矩阵计算公式是[ MAC(i,j) \frac{|\phi_i^T \phi_j|^2}{(\phi_i^T \phi_i)(\phi_j^T \phi_j)} ]对角元接近 1、非对角元接近 0 说明振型分离得干净。如果两个模态的 MAC 值超过 0.7基本上可以判定识别出了重复模态或振型混叠这时需要增加测点或重新选参考点。环境激励识别的阻尼比通常偏大原因是互相关函数尾部信噪比低ERA 拟合时会把噪声衰减也算进去。对策是把互相关函数截断到衰减至峰值 10% 左右的长度再喂给 ERA能修正约 0.3~0.5 个百分点的阻尼偏差。MATLAB 里直接取 R 序列前 1/3 段即可代价是频率分辨率略有下降。本文还有配套的精品资源点击获取