简介面向结构动力学与健康监测研究者的数据驱动随机子空间识别DDRSI算法实现包解决实际测量数据中结构模态参数自动提取问题适用于桥梁、建筑、航天器等设施的模态分析。压缩包共24个文件以13个MAT数据文件和11个M脚本为主M脚本覆盖稳定图绘制、模态参数提取、子空间计算等关键环节MAT文件则提供ASCE、Z24等benchmark数据及多组实测加速度输入便于按模块调用与替换数据。资源包27.33MB已有1116人学习/下载。读者可借助主程序串联整个流程结合多自由度连续结构等示例快速上手并利用自动识别模块直接获得频率、阻尼比与振型结果省去手动筛选稳定图的繁琐步骤适合具备一定MATLAB基础的结构工程与信号处理方向学习者。1. 为什么结构健康监测需要“自动”的模态识别在结构健康监测SHM里固有频率、阻尼比和振型是判断结构刚度退化最直接的证据。传统峰值拾取法和频域分解法在密集模态、强噪声场景下精度有限且依赖人工在频谱上圈选共振峰面对数万段在线监测数据人工成本不可接受。数据驱动的随机子空间识别DDRSI把加速度响应视为状态空间系统输出直接辨识一次SVD得到全阶候选模态再由稳定图自动筛选物理模态无需人工介入。下面基于SSI_DataDriven2代码包从Hankel矩阵构造开始把LQ分解、SVD截断、模态提取、稳定图自动聚类这条链路逐一拆解给出可直接运行的MATLAB实现和基于ASCE基准模型的调试参数。2. 随机子空间识别算法的理论基础Hankel矩阵与LQ投影2.1 从状态空间模型到Hankel矩阵的行空间关系DDRSI的理论起点是假设结构响应满足离散状态空间方程x_{k1} A x_k w_ky_k C x_k v_k其中 x_k 是 n 维不可观测状态y_k 是 m 维加速度响应w_k 是过程噪声v_k 是测量噪声。与传统时域AR模型不同随机子空间算法不直接拟合 y_k 的自回归系数而是把连续采样数组装成块Hankel矩阵利用“过去”与“未来”两个行空间之间的统计相关性来恢复系统的可观测性结构。为什么行空间这么重要结构在实际激励下的响应是确定性分量与外噪声的叠加确定性分量满足低阶状态空间关系噪声则没有这种关系。如果把未来响应投影到过去响应的行空间上与过去线性相关的那部分被保留下来不相关的随机分量被丢弃——这就是为什么随机子空间在强噪声环境下仍能稳定识别频率投影操作本身完成了第一轮信噪分离。2.2 blkhank.mHankel矩阵构造与块行数参数选择包内 blkhank.m 对应的功能实现如下function H blkhank(y, i) % 构造数据驱动随机子空间识别的块Hankel矩阵 % 输入: % y : m×N 加速度响应矩阵, m为测点通道数, N为每通道采样点数 % i : 块行数, 要求矩阵总行数 2*i*m 远小于 N % 输出: % H : 2*i*m × (N-2*i1) 的块Hankel矩阵 [m, N] size(y); cols N - 2*i 1; % 可用的时间窗口数 H zeros(2*i*m, cols); for k 1:2*i H((k-1)*m1 : k*m, :) y(:, k : kcols-1); end end循环里每推移一个块行截取的时窗就右移一个采样点最终形成 2i 个块行。前 i 个块行y 的第 1 到 i 段在算法里称为“过去”后 i 个块行y 的第 i1 到 2i 段称为“未来”。参数 i 直接决定矩阵规模可识别的系统阶数上限由 i*m 决定实际使用中 i 取太大时Hankel 矩阵接近病态SVD 的奇异值差距被抹平稳定图上会出现大量伪模态。我一般按两个经验值选 i一是保证 cols N-2i1 ≥ 500让后续 SVD 统计充分二是 2im ≤ N/3防止矩阵行数过大消耗内存。对 N120000、8 通道的加速度记录i20~40 是常用范围。2.3 LQ分解与斜投影SSI_CVA_Alg3.m的数学内核把块Hankel矩阵上下分块后记 Yp 为过去行块Yf 为未来行块。对 [Yp; Yf] 做 LQ 分解[Yp; Yf] L * Q其中 L 是下三角Q 是正交矩阵。分块后L [L11 0 ; L21 L22]Q [Q1 ; Q2]未来行块 Yf 在过去行空间 Yp 上的斜投影定义为 O Yf / Yp其 LQ 表达式就是 L21 与 Q1 的乘积。投影结果中噪声贡献被压缩到 L22*Q2 的残差项里之后对投影矩阵 O 做 SVD即可得到可观测矩阵 Gamma 的估计。SSI_CVA_Alg3.m 是包内的一组变体CVACanonical Variate Analysis在投影前后对过去和未来行做加权归一化。加权矩阵的选择会影响奇异值分布——经典数据驱动SSI直接对 O 做SVDCVA 把 L11 和 L22 的逆参与进来把投影相对能量差异拉平密集模态的分离度更好但对噪声也更敏感。我的建议是测点数少、信噪比高时优先CVA测点数据量大、噪声水平高时用不加权的投影会更稳。实际计算代码与包内 SSI_CVA_Alg3.m 思路一致function [A, C] DDRSI_Alg(y, i, n) % 数据驱动随机子空间核心计算 % y : m×N 加速度响应 % i : Hankel块行数 % n : 截断阶数 [m, N] size(y); H blkhank(y, i); p i*m; Yp H(1:p, :); Yf H(p1:2*p, :); % LQ分解: 对转置矩阵做经济型QR [Qr, Rr] qr([Yp; Yf], 0); L Rr; % 下三角 Q Qr; % 正交矩阵 L21 L(p1:2*p, 1:p); Q1 Q(1:p, :); % 斜投影 SVD O L21 * Q1; [U, S, ~] svd(O, econ); % 截断并恢复可观测矩阵 Gamma U(:, 1:n) * sqrt(S(1:n, 1:n)); % 位移不变性恢复系统矩阵 C Gamma(1:m, :); A Gamma(1:(i-1)*m, :) \ Gamma(m1:i*m, :); end逻辑说明qr 作用于转置矩阵得到 L 和 Q 后利用 L21 与 Q1 直接计算投影避免显式计算协方差矩阵这是数据驱动方法与协方差驱动方法最主要的区别。A 的恢复利用 Gamma 的位移不变性去掉最上面的 m 行与去掉最下面的 m 行之间满足线性映射关系最小二乘一次解出离散系统矩阵 A。参数说明n 是截断的奇异值个数对应系统阶数结构模态以共轭复极点成对出现所以 n 通常取偶数。包里的 solvric.m 和 predict.m 属于另一条实现路径——用稳态卡尔曼滤波器重构状态后再估计系统矩阵如果走本文这条投影SVD路线这两个脚本不需要介入。注意solvric.m 和 predict.m 与投影SVD路线不能混用。前者先求解Riccati方程获得卡尔曼增益再通过预测残差修正状态后者直接从投影分解拿状态。两者输出状态的含义不同混用会得到错误的A矩阵。3. 从SVD截断到模态参数系统矩阵恢复与特征值分解3.1 离散特征值为什么必须做对数逆映射DDRSI_Alg 输出的 A 是离散状态矩阵它的特征值 lambda_d 与连续时间特征值 lambda_c 满足 lambda_d exp(lambda_c * dt)。直接用离散特征值报频率会得到一个随采样率漂移的数值工程上要的是物理频率和阻尼比所以必须做一次对数逆映射lambda_c log(lambda_d) / dt这个映射对低阻尼结构的高频段尤其重要。采样率 fs200Hz、模态频率 30Hz 时离散特征值接近单位圆的高频边界一阶泰勒展开的误差会被放大只有严格对数映射才能保证阻尼比识别误差控制在可接受范围。注意低阻尼结构阻尼比小于 1%里对数映射的一阶近似误差虽然对频率影响小却会通过实部直接污染阻尼比。阻尼比数值本来就小相对误差放大两三个数量级都不奇怪所以不要省略这个映射步骤。3.2 GetMode_SSI.m频率、阻尼比与振型一次取回直接给出 GetMode_SSI.m 的等价实现function [fn, zeta, Phi] GetMode_SSI(A, C, dt) % 从离散状态空间矩阵提取模态参数 % A : n×n 离散系统矩阵 % C : m×n 输出矩阵 % dt: 采样间隔, 单位秒 [Psi, D] eig(A); lambda_d diag(D); lambda_c log(lambda_d) / dt; % 离散-连续映射 fn abs(lambda_c) / (2*pi); % 固有频率, Hz zeta -real(lambda_c) ./ abs(lambda_c); % 阻尼比 Phi C * Psi; % 测点上的振型分量 % 按频率升序排序输出 [fn, idx] sort(fn); zeta zeta(idx); Phi Phi(:, idx); end逻辑说明eig(A) 返回特征向量矩阵 Psi 和对角特征值 D。特征向量的物理含义是该阶模态在状态空间里的“形变方向”与输出矩阵 C 相乘后投影为 m 个测点上的振型分量。阻尼比是连续特征值实部与模长的负比值取绝对值是为了兼容可能出现的正实部噪声极点后续通过稳定图过滤。参数说明dt 一定取实际采样间隔不是采样频率的倒数这个顺序搞反会让全部频率和阻尼比错乱。输出按频率升序排序方便与有限元计算结果或峰值拾取法结果逐阶对照。3.3 阶数扫描策略为什么步长固定为2系统真实阶数未知工程上通过扫描 n 并观察模态参数随 n 的变化来筛选真实模态扫描参数范围说明起始阶数2最低一阶复模态对步长2每次加入一对共轭极点最大阶数不超过 i*m且保持偶数受Hankel块行数约束相邻阶判据频率变化≤1%阻尼变化≤5%MAC≥0.95稳定点定义为什么扫描步长固定为 2因为结构模态在状态空间里总是以共轭复极点对的形式出现每加入一对极点系统阶数增加 2。步长取 2 可以保证扫描轨迹与极点对结构对齐。另一个容易被忽视的约束是 i 的下限。当 i 小于真实系统阶数的一半时Hankel 矩阵装不下完整的可观性结构最大扫描阶数不足稳定图会出现“断层”——部分真实模态在高阶区域完全缺失。如果发现稳定图高阶区空白一段后又有新的极点出现先检查 i而不是去调容差参数。每一组 n 计算出的候选模态按“频率-阻尼”坐标画出来就是下一章要用到的稳定图。4. 稳定图与模态自动识别容差配置与SBMode_SSI.m聚类4.1 稳定点判定频率、阻尼、MAC三层判据衔接第二章的DDRSI_Alg输出稳定点判定的核心逻辑如下% 稳定点三层判据对应 SSI_StablizationDiagram.m 核心逻辑 function [f_stable, z_stable, flag] judge_stable(... f_cur, z_cur, Phi_cur, f_ref, z_ref, Phi_ref, tol) f_stable []; z_stable []; for k 1:length(f_cur) for j 1:length(f_ref) % 第一层: 频率相对误差 if abs(f_cur(k) - f_ref(j)) / f_ref(j) tol.freq continue; end % 第二层: 阻尼比绝对误差 if abs(z_cur(k) - z_ref(j)) tol.damp continue; end % 第三层: 振型MAC MAC abs(Phi_cur(:,k) * Phi_ref(:,j))^2 / ... ( (Phi_cur(:,k) * Phi_cur(:,k)) * (Phi_ref(:,j) * Phi_ref(:,j)) ); if MAC tol.MAC continue; end f_stable(end1) f_cur(k); z_stable(end1) z_cur(k); break; end end flag ~isempty(f_stable); end逻辑说明三层判断顺序不可颠倒。频率判断计算量最低先筛掉绝大多数无关极点阻尼对真实模态非常敏感识别方差天然比频率大容差要放宽MAC 放在最后即便频率和阻尼都接近的噪声极点只要振型形状不对就被刷掉。一个候选极点只允许匹配一个参考极点避免一次匹配引发连锁错误。参数说明这里的 f_ref、z_ref、Phi_ref 来自上一阶次的稳定极点集。实际工程中更推荐的写法是“先找最近频率再验证阻尼和MAC”而不是按固定顺序遍历所有参考极点后者的计算开销在阶次达到 60 时会有明显增长。4.2 容差参数怎么调一张表说清楚参数名默认值什么时候收紧/放宽TolFreq 频率相对误差1%模态间隔小于3%时收紧到0.5%强噪声下放宽到2%TolDamp 阻尼绝对误差0.05钢混结构阻尼本底0.01~0.030.05足够强风响应放宽到0.08TolMAC 振型一致性0.98测点密、通道多时0.98测点稀疏时放宽到0.90~0.95连续稳定阶数3Z24桥实测数据连续2阶即可接受阻尼比容差是最容易踩坑的参数。实验室环境里模态阻尼比识别方差很小0.05 的绝对容差绰绰有余但工程实测时环境激励的非平稳性会让同一模态在不同数据段的阻尼识别结果在 0.01 到 0.04 之间波动“稳定点”数量骤减。这时候先把 TolDamp 放宽到 0.08再看稳定图形态是否改善比反复调频率容差更有效。4.3 从稳定点到自动化输出滑窗聚类实现稳定图上满足判据的极点仍然是离散点需要再次聚类才能输出“一阶一值”的最终模态参数。SBMode_SSI.m 的做法可以概括为滑窗计数% 滑窗聚类简化自 SBMode_SSI.m function [f_final, z_final] cluster_stable(f_all, z_all, width) % f_all: 所有稳定点频率 % z_all: 对应阻尼比 % width: 聚类窗口宽度, 一般取最小频率间隔的一半 [f_sorted, idx] sort(f_all); z_sorted z_all(idx); f_final []; z_final []; while ~isempty(f_sorted) peak f_sorted(1); % 窗口内首个频率 in_win abs(f_sorted - peak) width; % 窗口内所有点 f_final(end1) median(f_sorted(in_win)); % 中位数代表该簇频率 z_final(end1) median(z_sorted(in_win)); % 中位数代表该簇阻尼 f_sorted(in_win) []; z_sorted(in_win) []; end end逻辑说明在自动模态识别中聚类窗口宽度的选择比判据容差更关键。取值为最小相邻真实模态频率间隔的一半是最安全的起点。对ASCE基准数据前三阶频率间隔在 1~3Hz窗口取 0.5Hz 能得到三个清晰的簇对密集模态结构窗口要缩小到 0.1Hz否则两个真实模态会被合并成一个伪模态。参数说明用中位数而不是均值作为簇代表原因是稳定点集合里混有少数未被容差完全过滤的离群噪声点中位数对离群值不敏感。当稳定点总数少于阶数的三倍时聚类结果不可信此时应放宽 TolDamp 或降低连续稳定阶数要求先把稳定点数提上来。5. ASCE基准模型与Z24实测数据的验证噪声、参数与三个坑5.1 同一组参数跑三条噪声数据集包内自带 ASCE(d0.02)_NoNoise.mat、ASCE(d0.02)_Noise5.mat、ASCE(d0.02)_Noise20.mat三组数据来自同一基准结构的同一工况区别只在信噪比。验证流程固定同一组 i20、n30对比前三阶频率% 统一参数对比三条数据集 files {ASCE(d0.02)_NoNoise.mat, ... ASCE(d0.02)_Noise5.mat, ... ASCE(d0.02)_Noise20.mat}; for k 1:length(files) S load(files{k}); fn_list fieldnames(S); y S.(fn_list{1}); % 取mat内第一个变量作为响应矩阵 if size(y, 1) size(y, 2) % 确保是 通道×采样点 y y; end [A, C] DDRSI_Alg(y, 20, 30); [fn, ~, ~] GetMode_SSI(A, C, 1/200); fprintf(%-28s f1%.4f f2%.4f f3%.4f\n, files{k}, fn(1), fn(2), fn(3)); end逻辑说明把阶数固定在 n30 而不是用稳定图扫描是为了让三条数据在完全相同的条件下对比。无噪声数据的前三阶频率与20%噪声数据的前三阶频率差距通常不超过 2%如果偏差超过这个幅度问题大概率出在块行数 i 或 Hankel 矩阵列数不足而不是算法本身。fs200Hz 对应的 dt1/200 来自包内数据文件的采样率约定。5.2 三个常见坑坑一测点通道数太少时 CVA 加权会失效。通道数 m 只有 2~3 个时SSI_CVA_Alg3.m 的加权矩阵接近奇异投影结果受最小奇异值支配稳定图出现大量高水平谱线。此时应该退回不加权投影或者把 i 至少加到 30 以上补偿信息量。坑二Z24桥实测数据要单独调 TolDamp。工程实测数据与实验室基准模型的差异集中体现在阻尼上实测阻尼比经常在 0.8%~1.5% 之间波动且随环境温度变化。把 TolDamp 从 0.05 放宽到 0.08 才能让稳定图连续频率容差反而要收紧到 0.5%因为实测频率对采样率和温度变化的敏感性更高。坑三连续稳定阶数与“模态真实存在”不能画等号。噪声白化特性好时噪声极点也可能形成连续 3 阶甚至 5 阶的稳定点。判断最终模态时把自动识别结果与 SVD 奇异值谱做一次交叉验证真实模态对应的奇异值通常在谱上形成明显拐点噪声模态的奇异值沿平滑曲线衰减。这一条验证步骤比任何自动聚类参数都更能防止误报模态。本文还有配套的精品资源点击获取 数据驱动随机子空间识别实现结构自动模态识别
简介面向结构动力学与健康监测研究者的数据驱动随机子空间识别DDRSI算法实现包解决实际测量数据中结构模态参数自动提取问题适用于桥梁、建筑、航天器等设施的模态分析。压缩包共24个文件以13个MAT数据文件和11个M脚本为主M脚本覆盖稳定图绘制、模态参数提取、子空间计算等关键环节MAT文件则提供ASCE、Z24等benchmark数据及多组实测加速度输入便于按模块调用与替换数据。资源包27.33MB已有1116人学习/下载。读者可借助主程序串联整个流程结合多自由度连续结构等示例快速上手并利用自动识别模块直接获得频率、阻尼比与振型结果省去手动筛选稳定图的繁琐步骤适合具备一定MATLAB基础的结构工程与信号处理方向学习者。1. 为什么结构健康监测需要“自动”的模态识别在结构健康监测SHM里固有频率、阻尼比和振型是判断结构刚度退化最直接的证据。传统峰值拾取法和频域分解法在密集模态、强噪声场景下精度有限且依赖人工在频谱上圈选共振峰面对数万段在线监测数据人工成本不可接受。数据驱动的随机子空间识别DDRSI把加速度响应视为状态空间系统输出直接辨识一次SVD得到全阶候选模态再由稳定图自动筛选物理模态无需人工介入。下面基于SSI_DataDriven2代码包从Hankel矩阵构造开始把LQ分解、SVD截断、模态提取、稳定图自动聚类这条链路逐一拆解给出可直接运行的MATLAB实现和基于ASCE基准模型的调试参数。2. 随机子空间识别算法的理论基础Hankel矩阵与LQ投影2.1 从状态空间模型到Hankel矩阵的行空间关系DDRSI的理论起点是假设结构响应满足离散状态空间方程x_{k1} A x_k w_ky_k C x_k v_k其中 x_k 是 n 维不可观测状态y_k 是 m 维加速度响应w_k 是过程噪声v_k 是测量噪声。与传统时域AR模型不同随机子空间算法不直接拟合 y_k 的自回归系数而是把连续采样数组装成块Hankel矩阵利用“过去”与“未来”两个行空间之间的统计相关性来恢复系统的可观测性结构。为什么行空间这么重要结构在实际激励下的响应是确定性分量与外噪声的叠加确定性分量满足低阶状态空间关系噪声则没有这种关系。如果把未来响应投影到过去响应的行空间上与过去线性相关的那部分被保留下来不相关的随机分量被丢弃——这就是为什么随机子空间在强噪声环境下仍能稳定识别频率投影操作本身完成了第一轮信噪分离。2.2 blkhank.mHankel矩阵构造与块行数参数选择包内 blkhank.m 对应的功能实现如下function H blkhank(y, i) % 构造数据驱动随机子空间识别的块Hankel矩阵 % 输入: % y : m×N 加速度响应矩阵, m为测点通道数, N为每通道采样点数 % i : 块行数, 要求矩阵总行数 2*i*m 远小于 N % 输出: % H : 2*i*m × (N-2*i1) 的块Hankel矩阵 [m, N] size(y); cols N - 2*i 1; % 可用的时间窗口数 H zeros(2*i*m, cols); for k 1:2*i H((k-1)*m1 : k*m, :) y(:, k : kcols-1); end end循环里每推移一个块行截取的时窗就右移一个采样点最终形成 2i 个块行。前 i 个块行y 的第 1 到 i 段在算法里称为“过去”后 i 个块行y 的第 i1 到 2i 段称为“未来”。参数 i 直接决定矩阵规模可识别的系统阶数上限由 i*m 决定实际使用中 i 取太大时Hankel 矩阵接近病态SVD 的奇异值差距被抹平稳定图上会出现大量伪模态。我一般按两个经验值选 i一是保证 cols N-2i1 ≥ 500让后续 SVD 统计充分二是 2im ≤ N/3防止矩阵行数过大消耗内存。对 N120000、8 通道的加速度记录i20~40 是常用范围。2.3 LQ分解与斜投影SSI_CVA_Alg3.m的数学内核把块Hankel矩阵上下分块后记 Yp 为过去行块Yf 为未来行块。对 [Yp; Yf] 做 LQ 分解[Yp; Yf] L * Q其中 L 是下三角Q 是正交矩阵。分块后L [L11 0 ; L21 L22]Q [Q1 ; Q2]未来行块 Yf 在过去行空间 Yp 上的斜投影定义为 O Yf / Yp其 LQ 表达式就是 L21 与 Q1 的乘积。投影结果中噪声贡献被压缩到 L22*Q2 的残差项里之后对投影矩阵 O 做 SVD即可得到可观测矩阵 Gamma 的估计。SSI_CVA_Alg3.m 是包内的一组变体CVACanonical Variate Analysis在投影前后对过去和未来行做加权归一化。加权矩阵的选择会影响奇异值分布——经典数据驱动SSI直接对 O 做SVDCVA 把 L11 和 L22 的逆参与进来把投影相对能量差异拉平密集模态的分离度更好但对噪声也更敏感。我的建议是测点数少、信噪比高时优先CVA测点数据量大、噪声水平高时用不加权的投影会更稳。实际计算代码与包内 SSI_CVA_Alg3.m 思路一致function [A, C] DDRSI_Alg(y, i, n) % 数据驱动随机子空间核心计算 % y : m×N 加速度响应 % i : Hankel块行数 % n : 截断阶数 [m, N] size(y); H blkhank(y, i); p i*m; Yp H(1:p, :); Yf H(p1:2*p, :); % LQ分解: 对转置矩阵做经济型QR [Qr, Rr] qr([Yp; Yf], 0); L Rr; % 下三角 Q Qr; % 正交矩阵 L21 L(p1:2*p, 1:p); Q1 Q(1:p, :); % 斜投影 SVD O L21 * Q1; [U, S, ~] svd(O, econ); % 截断并恢复可观测矩阵 Gamma U(:, 1:n) * sqrt(S(1:n, 1:n)); % 位移不变性恢复系统矩阵 C Gamma(1:m, :); A Gamma(1:(i-1)*m, :) \ Gamma(m1:i*m, :); end逻辑说明qr 作用于转置矩阵得到 L 和 Q 后利用 L21 与 Q1 直接计算投影避免显式计算协方差矩阵这是数据驱动方法与协方差驱动方法最主要的区别。A 的恢复利用 Gamma 的位移不变性去掉最上面的 m 行与去掉最下面的 m 行之间满足线性映射关系最小二乘一次解出离散系统矩阵 A。参数说明n 是截断的奇异值个数对应系统阶数结构模态以共轭复极点成对出现所以 n 通常取偶数。包里的 solvric.m 和 predict.m 属于另一条实现路径——用稳态卡尔曼滤波器重构状态后再估计系统矩阵如果走本文这条投影SVD路线这两个脚本不需要介入。注意solvric.m 和 predict.m 与投影SVD路线不能混用。前者先求解Riccati方程获得卡尔曼增益再通过预测残差修正状态后者直接从投影分解拿状态。两者输出状态的含义不同混用会得到错误的A矩阵。3. 从SVD截断到模态参数系统矩阵恢复与特征值分解3.1 离散特征值为什么必须做对数逆映射DDRSI_Alg 输出的 A 是离散状态矩阵它的特征值 lambda_d 与连续时间特征值 lambda_c 满足 lambda_d exp(lambda_c * dt)。直接用离散特征值报频率会得到一个随采样率漂移的数值工程上要的是物理频率和阻尼比所以必须做一次对数逆映射lambda_c log(lambda_d) / dt这个映射对低阻尼结构的高频段尤其重要。采样率 fs200Hz、模态频率 30Hz 时离散特征值接近单位圆的高频边界一阶泰勒展开的误差会被放大只有严格对数映射才能保证阻尼比识别误差控制在可接受范围。注意低阻尼结构阻尼比小于 1%里对数映射的一阶近似误差虽然对频率影响小却会通过实部直接污染阻尼比。阻尼比数值本来就小相对误差放大两三个数量级都不奇怪所以不要省略这个映射步骤。3.2 GetMode_SSI.m频率、阻尼比与振型一次取回直接给出 GetMode_SSI.m 的等价实现function [fn, zeta, Phi] GetMode_SSI(A, C, dt) % 从离散状态空间矩阵提取模态参数 % A : n×n 离散系统矩阵 % C : m×n 输出矩阵 % dt: 采样间隔, 单位秒 [Psi, D] eig(A); lambda_d diag(D); lambda_c log(lambda_d) / dt; % 离散-连续映射 fn abs(lambda_c) / (2*pi); % 固有频率, Hz zeta -real(lambda_c) ./ abs(lambda_c); % 阻尼比 Phi C * Psi; % 测点上的振型分量 % 按频率升序排序输出 [fn, idx] sort(fn); zeta zeta(idx); Phi Phi(:, idx); end逻辑说明eig(A) 返回特征向量矩阵 Psi 和对角特征值 D。特征向量的物理含义是该阶模态在状态空间里的“形变方向”与输出矩阵 C 相乘后投影为 m 个测点上的振型分量。阻尼比是连续特征值实部与模长的负比值取绝对值是为了兼容可能出现的正实部噪声极点后续通过稳定图过滤。参数说明dt 一定取实际采样间隔不是采样频率的倒数这个顺序搞反会让全部频率和阻尼比错乱。输出按频率升序排序方便与有限元计算结果或峰值拾取法结果逐阶对照。3.3 阶数扫描策略为什么步长固定为2系统真实阶数未知工程上通过扫描 n 并观察模态参数随 n 的变化来筛选真实模态扫描参数范围说明起始阶数2最低一阶复模态对步长2每次加入一对共轭极点最大阶数不超过 i*m且保持偶数受Hankel块行数约束相邻阶判据频率变化≤1%阻尼变化≤5%MAC≥0.95稳定点定义为什么扫描步长固定为 2因为结构模态在状态空间里总是以共轭复极点对的形式出现每加入一对极点系统阶数增加 2。步长取 2 可以保证扫描轨迹与极点对结构对齐。另一个容易被忽视的约束是 i 的下限。当 i 小于真实系统阶数的一半时Hankel 矩阵装不下完整的可观性结构最大扫描阶数不足稳定图会出现“断层”——部分真实模态在高阶区域完全缺失。如果发现稳定图高阶区空白一段后又有新的极点出现先检查 i而不是去调容差参数。每一组 n 计算出的候选模态按“频率-阻尼”坐标画出来就是下一章要用到的稳定图。4. 稳定图与模态自动识别容差配置与SBMode_SSI.m聚类4.1 稳定点判定频率、阻尼、MAC三层判据衔接第二章的DDRSI_Alg输出稳定点判定的核心逻辑如下% 稳定点三层判据对应 SSI_StablizationDiagram.m 核心逻辑 function [f_stable, z_stable, flag] judge_stable(... f_cur, z_cur, Phi_cur, f_ref, z_ref, Phi_ref, tol) f_stable []; z_stable []; for k 1:length(f_cur) for j 1:length(f_ref) % 第一层: 频率相对误差 if abs(f_cur(k) - f_ref(j)) / f_ref(j) tol.freq continue; end % 第二层: 阻尼比绝对误差 if abs(z_cur(k) - z_ref(j)) tol.damp continue; end % 第三层: 振型MAC MAC abs(Phi_cur(:,k) * Phi_ref(:,j))^2 / ... ( (Phi_cur(:,k) * Phi_cur(:,k)) * (Phi_ref(:,j) * Phi_ref(:,j)) ); if MAC tol.MAC continue; end f_stable(end1) f_cur(k); z_stable(end1) z_cur(k); break; end end flag ~isempty(f_stable); end逻辑说明三层判断顺序不可颠倒。频率判断计算量最低先筛掉绝大多数无关极点阻尼对真实模态非常敏感识别方差天然比频率大容差要放宽MAC 放在最后即便频率和阻尼都接近的噪声极点只要振型形状不对就被刷掉。一个候选极点只允许匹配一个参考极点避免一次匹配引发连锁错误。参数说明这里的 f_ref、z_ref、Phi_ref 来自上一阶次的稳定极点集。实际工程中更推荐的写法是“先找最近频率再验证阻尼和MAC”而不是按固定顺序遍历所有参考极点后者的计算开销在阶次达到 60 时会有明显增长。4.2 容差参数怎么调一张表说清楚参数名默认值什么时候收紧/放宽TolFreq 频率相对误差1%模态间隔小于3%时收紧到0.5%强噪声下放宽到2%TolDamp 阻尼绝对误差0.05钢混结构阻尼本底0.01~0.030.05足够强风响应放宽到0.08TolMAC 振型一致性0.98测点密、通道多时0.98测点稀疏时放宽到0.90~0.95连续稳定阶数3Z24桥实测数据连续2阶即可接受阻尼比容差是最容易踩坑的参数。实验室环境里模态阻尼比识别方差很小0.05 的绝对容差绰绰有余但工程实测时环境激励的非平稳性会让同一模态在不同数据段的阻尼识别结果在 0.01 到 0.04 之间波动“稳定点”数量骤减。这时候先把 TolDamp 放宽到 0.08再看稳定图形态是否改善比反复调频率容差更有效。4.3 从稳定点到自动化输出滑窗聚类实现稳定图上满足判据的极点仍然是离散点需要再次聚类才能输出“一阶一值”的最终模态参数。SBMode_SSI.m 的做法可以概括为滑窗计数% 滑窗聚类简化自 SBMode_SSI.m function [f_final, z_final] cluster_stable(f_all, z_all, width) % f_all: 所有稳定点频率 % z_all: 对应阻尼比 % width: 聚类窗口宽度, 一般取最小频率间隔的一半 [f_sorted, idx] sort(f_all); z_sorted z_all(idx); f_final []; z_final []; while ~isempty(f_sorted) peak f_sorted(1); % 窗口内首个频率 in_win abs(f_sorted - peak) width; % 窗口内所有点 f_final(end1) median(f_sorted(in_win)); % 中位数代表该簇频率 z_final(end1) median(z_sorted(in_win)); % 中位数代表该簇阻尼 f_sorted(in_win) []; z_sorted(in_win) []; end end逻辑说明在自动模态识别中聚类窗口宽度的选择比判据容差更关键。取值为最小相邻真实模态频率间隔的一半是最安全的起点。对ASCE基准数据前三阶频率间隔在 1~3Hz窗口取 0.5Hz 能得到三个清晰的簇对密集模态结构窗口要缩小到 0.1Hz否则两个真实模态会被合并成一个伪模态。参数说明用中位数而不是均值作为簇代表原因是稳定点集合里混有少数未被容差完全过滤的离群噪声点中位数对离群值不敏感。当稳定点总数少于阶数的三倍时聚类结果不可信此时应放宽 TolDamp 或降低连续稳定阶数要求先把稳定点数提上来。5. ASCE基准模型与Z24实测数据的验证噪声、参数与三个坑5.1 同一组参数跑三条噪声数据集包内自带 ASCE(d0.02)_NoNoise.mat、ASCE(d0.02)_Noise5.mat、ASCE(d0.02)_Noise20.mat三组数据来自同一基准结构的同一工况区别只在信噪比。验证流程固定同一组 i20、n30对比前三阶频率% 统一参数对比三条数据集 files {ASCE(d0.02)_NoNoise.mat, ... ASCE(d0.02)_Noise5.mat, ... ASCE(d0.02)_Noise20.mat}; for k 1:length(files) S load(files{k}); fn_list fieldnames(S); y S.(fn_list{1}); % 取mat内第一个变量作为响应矩阵 if size(y, 1) size(y, 2) % 确保是 通道×采样点 y y; end [A, C] DDRSI_Alg(y, 20, 30); [fn, ~, ~] GetMode_SSI(A, C, 1/200); fprintf(%-28s f1%.4f f2%.4f f3%.4f\n, files{k}, fn(1), fn(2), fn(3)); end逻辑说明把阶数固定在 n30 而不是用稳定图扫描是为了让三条数据在完全相同的条件下对比。无噪声数据的前三阶频率与20%噪声数据的前三阶频率差距通常不超过 2%如果偏差超过这个幅度问题大概率出在块行数 i 或 Hankel 矩阵列数不足而不是算法本身。fs200Hz 对应的 dt1/200 来自包内数据文件的采样率约定。5.2 三个常见坑坑一测点通道数太少时 CVA 加权会失效。通道数 m 只有 2~3 个时SSI_CVA_Alg3.m 的加权矩阵接近奇异投影结果受最小奇异值支配稳定图出现大量高水平谱线。此时应该退回不加权投影或者把 i 至少加到 30 以上补偿信息量。坑二Z24桥实测数据要单独调 TolDamp。工程实测数据与实验室基准模型的差异集中体现在阻尼上实测阻尼比经常在 0.8%~1.5% 之间波动且随环境温度变化。把 TolDamp 从 0.05 放宽到 0.08 才能让稳定图连续频率容差反而要收紧到 0.5%因为实测频率对采样率和温度变化的敏感性更高。坑三连续稳定阶数与“模态真实存在”不能画等号。噪声白化特性好时噪声极点也可能形成连续 3 阶甚至 5 阶的稳定点。判断最终模态时把自动识别结果与 SVD 奇异值谱做一次交叉验证真实模态对应的奇异值通常在谱上形成明显拐点噪声模态的奇异值沿平滑曲线衰减。这一条验证步骤比任何自动聚类参数都更能防止误报模态。本文还有配套的精品资源点击获取