1. SSI-COV方法值不值得学原理背景与选型判断先说结论如果你手里有一批结构在环境激励下的加速度响应数据想从中拿到模态频率、振型和阻尼比SSI-COV协方差驱动的随机子空间方法是目前最值得优先尝试的时域方法之一。我最近做一个结构健康监测方向的仿真验证项目对比了峰值拾取法、频域分解法和随机子空间法之后最后还是把SSI-COV作为主算法因为它在阻尼比识别和密集模态区分上的表现确实比频域那套要稳。1.1 为什么只靠响应就能识别模态传统实验模态分析你得有力锤或激振器测到输入力再算频响函数然后从频响函数里拟合模态参数。这个流程实验室里很好用但换到真实工程就尴尬了——一座跨江大桥、一栋超高层建筑你没法用激振器施加可控激励只能靠风、地脉动、车辆通行这类环境激励。环境激励的特点是激励力未知你手里只有传感器测到的响应数据。所以必须走仅输出的模态识别路线。仅输出识别的基本前提是环境激励可以近似看作随机白噪声或宽带平稳激励。在这个假设下响应的自相关函数和互相关函数里其实已经包含了系统的全部动力学信息。SSI-COV的核心思路就是先通过响应数据构造协方差序列再把这些协方差拼成Toeplitz矩阵最后对这个矩阵做奇异值分解把系统的状态矩阵从数据里投影出来。你不需要知道激励力怎么作用只需要保证激励的频带覆盖你关心的模态范围。1.2 随机子空间家族里SSI-COV的位置随机子空间方法内部还分两个方向SSI-COV协方差驱动和SSI-DATA数据驱动。SSI-DATA是从原始数据直接做LQ分解数值上更稳定但实现复杂度高矩阵维度一上去内存和计算量都不小。SSI-COV是先算协方差、再对Toeplitz矩阵做SVD思路更直观代码也更短非常适合自己动手实现和验证算法。代价是协方差估计这一步对数据长度有一定要求数据太短会导致协方差估计误差变大识别出来的阻尼比容易出现偏差。实际项目里怎么选如果数据记录足够长我一般要求至少覆盖感兴趣最低频率对应周期的50倍以上SSI-COV完全够用而且调试方便。如果你只有一小段数据或者信噪比很差那可以考虑SSI-DATA。本文后面所有内容都围绕SSI-COV展开因为把协方差驱动的原理吃透了再去看SSI-DATA就很轻松。提示模态参数识别里的阻尼比是最难识别准的参数不管哪种方法都逃不过数据长度和噪声水平的制约。后面你会看到识别频率误差能做到0.5%以内阻尼比误差却经常到5%~10%这是方法的固有特性不是代码写错了。2. 仿真算例三自由度系统的响应数据怎么来做算法验证的第一件事就是构造一个答案已知的系统。实际测量数据永远没有标准答案你无法判断算法识别出来的是对是错。所以我先用一个三自由度弹簧-质量-阻尼系统做仿真把理论模态参数先算出来再生成模拟响应数据最后用SSI-COV去识别看识别结果和理论值差多少。2.1 系统模型与状态空间表达我选的模型是一个串联形式的剪切型结构类似三层框架的简化模型。三个质量块都取1000kg层间刚度取4e6N/m阻尼采用瑞利阻尼[ \mathbf{M}m\mathbf{I},\quad \mathbf{K}k\begin{bmatrix}2 -1 0\ -1 2 -1\ 0 -1 1\end{bmatrix},\quad \mathbf{C}0.5\mathbf{M}0.0005\mathbf{K} ]质量、刚度矩阵都是典型的三自由度形式。阻尼用瑞利阻尼是为了方便计算理论阻尼比因为比例阻尼条件下模态阻尼比有解析表达式[ \zeta_i\frac{a}{2\omega_i}\frac{b\omega_i}{2} ]这里 (a0.5)、(b0.0005)可以让三个模态的阻尼比落在1.5%~3%这个比较真实的范围内不至于太小导致数值仿真难以稳定也不至于太大偏离实际结构。把运动方程改写成状态空间形式状态向量取位移和速度m 1000; k 4e6; M m * eye(3); K k * [2 -1 0; -1 2 -1; 0 -1 1]; C_damp 0.5 * M 0.0005 * K; % 状态空间连续时间矩阵 Ac [zeros(3), eye(3); -M\K, -M\C_damp]; Bc [zeros(3); inv(M)]; Cc [-M\K, -M\C_damp]; % 输出加速度 Dc inv(M); fs 200; dt 1/fs; t 0:dt:600; N length(t); f_ext randn(3, N); % 三个自由度上独立的随机激励 sys_c ss(Ac, Bc, Cc, Dc); y lsim(sys_c, f_ext, t); % 3 x N 的加速度响应采用加速度作为输出是因为工程现场绝大多数传感器是加速度计。注意这里Dc不是零当激励是作用于质量块上的力时加速度输出会直接包含力项。很多教程为了方便把D设成零那其实是假设激励通过某种方式不直接影响测量和实际情况有偏差。2.2 理论模态参数用于后续验证的标准答案三自由度系统的理论模态参数可以直接求解广义特征值问题[Phi_th, Om2] eig(K, M); [omega_th, idx] sort(sqrt(diag(Om2))); Phi_th Phi_th(:, idx); f_th omega_th / (2*pi); zeta_th 0.5 ./ (2*omega_th) 0.0005 * omega_th / 2;这套代码算出来三个模态的理论频率大约是4.48Hz、12.55Hz、18.13Hz理论阻尼比分别是1.59%、2.29%、3.07%。振型是三维向量后面跟识别结果做MAC对比时用。有一点要提醒这里的振型是位移振型而后面SSI-COV识别用的是加速度响应。加速度振型和位移振型之间差一个 (-\omega^2) 的缩放因子但因为每个模态都有各自的系数归一化之后两者是一致的做MAC相关性计算时不受影响。2.3 生成模拟加速度响应并加入噪声仿真激励用的是高斯白噪声目的是模拟环境脉动这类宽带随机激励。采样频率设200Hz数据长度600秒。加噪声时没有用固定的绝对噪声幅值而是按每个通道响应RMS的百分比加这样信噪比更可控noise_ratio 0.05; y y noise_ratio * std(y, 0, 2) .* randn(size(y));std(y,0,2)是求每个输出通道的标准差noise_ratio0.05表示噪声RMS是信号RMS的5%这个水平在仿真里已经不算干净了可以用来检验算法的抗噪能力。实际工程项目里现场数据的噪声水平经常比这更差所以下面还会单独分析不同噪声水平的影响。3. SSI-COV的Matlab实现算法拆解与代码逐段剖析这一章是整篇文章的核心。我会把SSI-COV算法的每个关键步骤拆开讲配上可以直接运行的Matlab代码。理解了每个矩阵在做什么你的算法调试能力会上一个台阶。3.1 Hankel矩阵、协方差矩阵与SVD截断SSI-COV虽然叫协方差驱动但工程实现时通常用Hankel矩阵的投影来计算效果等价、代码更简洁。先把响应数据排成Hankel矩阵分成过去和未来两块l size(y, 1); % 输出通道数这里是3 i_block 30; % 块行数 Ndata size(y, 2); ncol Ndata - 2*i_block 1; H zeros(2*i_block*l, ncol); for r 1:2*i_block H((r-1)*l1:r*l, :) y(:, r:rncol-1); end Yp H(1:i_block*l, :); % 过去块 Yf H(i_block*l1:2*i_block*l, :); % 未来块 T Yf * Yp / ncol; % l*i_block x l*i_block这个T矩阵就是算法要分解的对象。它的物理含义是把未来的响应数据往过去的响应数据方向上投影提取出由过去状态可预测的那部分结构信息剩下的部分属于新进入系统的随机激励被这一步过滤掉。接下来对该矩阵做奇异值分解[U, S, V] svd(T); n_order 8; % 先按略大于2倍物理模态数选取后面再讲稳定图定阶 O U(:, 1:n_order) * sqrt(S(1:n_order, 1:n_order));S矩阵的奇异值大小反映了各阶子空间在数据里的能量占比。系统只有3阶对应6个共轭极点理论上取n_order6就够。但实际数据有噪声、有随机误差SVD的奇异值不会干净地截断所以通常取大一点后面通过稳定图筛选真实模态。我这里的n_order8只是为了先把算法跑通。3.2 提取系统矩阵并换算频率、阻尼比、振型从截断后的可观测矩阵O可以提取输出矩阵C和系统矩阵A_d。关键关系是如果 (O[C;CA;CA^2;\cdots;CA^{i-1}])那么去掉前块和去掉后块的两部分满足 (\text{O_down}O_{up}\cdot A_d)。所以C_est O(1:l, :); A_d pinv(O(1:end-l, :)) * O(l1:end, :);A_d是离散时间状态矩阵要先做特征值分解[Psi, Lambda] eig(A_d); lambda_c log(diag(Lambda)) / dt;lambda_c就是连续时间系统的特征值是一对一对共轭复数。每一个共轭对对应一阶物理模态再进行参数换算fn_est abs(lambda_c) / (2*pi); zeta_est -real(lambda_c) ./ abs(lambda_c); Phi_est C_est * Psi;这里的物理逻辑是连续时间特征值 (\lambda-\zeta\omegaj\omega\sqrt{1-\zeta^2})它的模就是无阻尼固有圆频率实部的负值除以模就是阻尼比。振型向量是输出矩阵和特征向量的乘积 (C\Psi) 的各列按任意通道归一化后就是实际意义的模态振型。识别完成后把频率排序、剔除虚部为负的共轭重复项、只保留阻尼比在合理范围内的极点就得到最终结果。3.3 离散特征值到连续特征值的坑这一节必须单独拿出来说因为新手经常在这里翻车。离散状态矩阵的特征值 (\lambda_d) 和连续特征值 (\lambda_c) 的关系是[ \lambda_de^{\lambda_c\Delta t} ]所以反求连续特征值不能直接对矩阵开方而是要用矩阵对数lambda_c log(diag(Lambda)) / dt;Matlab里如果直接写logm(A_d)/dt也能算出结果但对角化情况下用log(diag(Lambda))更直观而且便于逐阶处理。麻烦的点在于复对数有分支问题如果采样频率不够高某些模态的离散特征值会落在单位圆的其它分支上导致换算出来的频率产生虚假偏移。我的经验是采样频率至少要覆盖感兴趣最高模态频率的5~10倍再低就容易出问题。此外识别结果里会出现成对的共轭极点频率完全相同这是正常现象。筛选取值时只保留其中一个即可千万别把共轭对当成两阶模态。提示阻尼比计算时如果lambda_c的实部为正说明识别出的极点不稳定对应的可能是噪声造成的伪模态可以直接丢弃。真实结构的阻尼比是正的但环境激励数据里偶尔会识别出负阻尼极点这在随机子空间方法里很常见。4. 结果验证频率、阻尼、振型到底准不准算法跑通了下一步就是看识别结果靠不靠谱。我用的是5%噪声、600秒数据的仿真算例随机种子固定后得到一组具体结果。下面从频率、阻尼比、振型三个维度逐一验证。4.1 模态参数对比三阶模态的识别结果和理论值对比如下模态理论频率 (Hz)识别频率 (Hz)频率误差理论阻尼比 (%)识别阻尼比 (%)阻尼误差1阶4.484.490.22%1.591.524.4%2阶12.5512.530.16%2.292.394.4%3阶18.1318.200.39%3.073.224.9%频率识别精度非常高三阶都控制在0.5%以内这符合SSI-COV的一贯表现。阻尼比的误差明显更大在4%~5%左右这也是所有仅输出法的通病。阻尼比的本质是能量耗散参数它隐含在响应衰减速度里对噪声、数据长度、频率分辨率都非常敏感。4.2 振型相关性MAC振型对比通常用模态置信准则MAC值。MAC的定义是在两个振型向量之间做一个相关性归一化[ MAC\frac{|\phi_1^H\phi_2|^2}{(\phi_1^H\phi_1)(\phi_2^H\phi_2)} ]数值越接近1说明两个振型越一致。Matlab里实现很简单function mac MAC(phi1, phi2) mac abs(phi1 * phi2)^2 / ((phi1 * phi1) * (phi2 * phi2)); end我用识别振型和理论振型逐阶做了MAC计算。三阶的MAC分别为0.9997、0.9991、0.9980说明振型识别得相当准。唯一要注意的是识别出的振型符号可能与理论振型相差180度这完全不影响MAC值因为模的平方会把符号消掉。提示做MAC时要注意向量数据类型。加速度响应的识别振量本质上跟位移振量差一个缩放但归一化后就不影响MAC了。实际工程中如果做多工况合并每个工况的振型也要先归一化再比较否则MAC计算就是白搭。4.3 噪声水平和数据长度的影响同一个算例我把噪声比从2%调高到10%再缩短数据长度观察表现。噪声在2%水平时频率误差小于0.1%阻尼误差约2%噪声升到10%时频率误差仍在1%以内但阻尼误差会扩大到10%~15%而且高阶模态可能出现伪极点需要靠稳定图手动筛选。数据长度的影响更值得注意。同样的5%噪声600秒数据识别出来的阻尼比误差约5%压缩到120秒后频率误差还在0.5%以内但阻尼比误差会飙到10%~20%。这说明如果目标是识别准确的阻尼比数据记录一定要足够长建议至少保证最低频率对应周期的50倍以上我个人习惯按100倍来取。5. 实战中绕不开的细节稳定图、预处理和定阶仿真环境里一切都干净利落真正做实测数据时各种问题才会暴露出来。这一章聊聊那些我在实际项目中反复踩过、花了不少时间才摸清楚的细节。5.1 模型阶次不是拍脑袋定的靠稳定图前面代码里我直接给了n_order8那是预设条件。实际数据里系统阶次完全未知而且你也不可能恰好知道结构有几阶模态在频带里。我的做法是把阶次从2循环到40对每个阶次都做一遍SSI-COV识别把所有识别出的极点画到一张频率-阶次图上这就是稳定图。稳定点的判断标准一般是相邻阶次之间频率变化小于1%、阻尼变化小于5%~10%、MAC大于0.95。满足这些条件的极点会在图上形成一条垂直的稳定线那条线对应的就是真实模态。手写稳定图代码不复杂核心逻辑是orderRange 2:2:40; stableFreq []; for n_order orderRange % 执行SSI-COV得到fn_est, zeta_est, Phi_est % 对每个候选极点和前面已确认的稳定极点比较 % 如果频率差1% 且 阻尼差5% 且 MAC0.95标记为稳定点 scatter(fn_est, n_order * ones(size(fn_est)), 20, fill); end千万不要直接选一个过高的阶次然后拿结果去讲物理故事那样会把噪声拟合得特别好但识别出来的模态毫无意义。5.2 数据预处理的正确姿势很多同学把原始数据拿过来就直接丢进SSI-COV结果识别出来一堆莫名其妙的小峰值。我踩过几次坑之后总结出一条基本流程第一去均值。加速度计现场记录数据经常有直流偏置不去均值的话协方差矩阵的第一行会有明显的趋势项污染识别结果在低频段会出现假模态。第二去趋势。数据本身如果有缓慢漂移先做多项式去趋势。我一般用一次或二次多项式趋势项对协方差估计的影响比均值更大。第三滤波。只在感兴趣的频带内保留信号比如你关心2~25Hz的模态就做2Hz高通和25Hz低通滤波。滤波能显著压低频带外噪声但要注意滤波器本身可能引入相位失真最好用零相位滤波函数如filtfilt不然识别出的阻尼比会被虚假地放大。第四重采样。原始数据采样率过高会让Hankel矩阵维度膨胀、计算量暴增。先低通再降采样到最高关注频率的5~10倍计算效率可以提升一个数量级。5.3 从仿真走向实测的三点提醒最后说三个很容易被忽略、却直接影响实测效果的点。传感器数量。SSI-COV是哪个通道都参与计算的。通道数太少你只能识别出测点能够分辨的振型高阶模态或者空间上相邻的模态会糊在一起。如果条件允许测点数量至少应该是你关心的模态阶数的一倍以上。激励的充分性。SSI-COV的白噪声激励假设在实际中永远是近似满足的。如果有风致涡激振动、设备转速引起的周期信号会在响应谱里叠加窄带峰值SSI会把这些也识别成模态。遇到这种情况要么换一段数据要么对这个频带做专门处理别硬套算法。振型归一化约定。做完振型识别在报告里说明你用的是哪种归一化是按第一个传感器位置归一还是按质量归一。不同归一化方式会直接影响后续损伤识别、模型修正的计算结果。我平时习惯统一按最大幅值归一跨工况对比时省很多纠结。我个人经验是SSI-COV这个算法调试一旦过了稳定图这道门槛后面就顺畅了。下一步如果想进阶可以试试在协方差的基础上加入加权矩阵或者把SSI-COV识别的状态矩阵拿去和有限元模型做相关性分析那又是另一个很有意思的方向了。