做结构振动测试的人十有八九遇到过这种尴尬测了一堆加速度响应数据激励力却测不到。风、车流、环境微振动一直在激励结构我们根本没有办法记录真实的输入。这时候想做模态参数识别传统的频响函数法直接失效只能用“仅输出”的方法。而在这一众方法里基于SSI-COV协方差驱动随机子空间识别的多自由度系统模态参数识别是我在Matlab里用得最顺手、也最愿意跟人推荐的一套。它不依赖人工激励不用迭代求解从响应数据的协方差统计量出发走一遍SVD和特征值分解就能把模态频率、振型和阻尼比全部提出来。这篇文章想把SSI-COV从原理到代码到工程注意事项完整讲透。内容主要面向做结构健康监测、振动测试、动力学仿真的研究生和工程师同时也适合刚接触模态识别、想在Matlab里快速落地一套可用流程的初学者。文中会包含可直接运行的Matlab核心代码、参数选取经验、稳定图实现思路以及一个三自由度系统的完整仿真验证算例。看完之后你可以直接拿这套东西去处理自己的实测数据也可以根据自己的问题场景做二次开发。1. 为什么单选SSI-COV环境激励下的识别困局1.1 没激励信号时传统方法为什么玩不转经典模态分析里最常用的方法是同时采集输入力信号和结构响应信号然后算频响函数FRF再从幅频曲线的峰值去对频率和阻尼。这个方法在实验室里非常成熟激振器一装、力锤一敲FRF干净漂亮参数识别精度很高。但到了现场问题就来了风荷载多大车流激励怎么测环境微振动压根没有一个明确的激励源。你既不可能在桥墩上装力传感器也不可能把整个结构停下来单独做一次受迫振动试验。这种情况下基于输入输出的方法全部失效。于是业内转向“仅输出”模态识别我不需要知道激励是什么只假设它是平稳随机信号从结构响应的统计规律里把模态信息反推出来。常用的仅输出方法有几类。频域分解法FDD思路简单把响应功率谱密度矩阵做SVD奇异值曲线的峰值对应模态频率但它给出的阻尼比精度一般只能作为初步参考。随机减量法RDT或NExT技术需要先从响应里提取自由衰减响应再做时域拟合可是提取过程中很依赖经验参数稍不注意自由衰减信号就变形了。Eigensystem Realization AlgorithmERA本身是个好方法但它需要脉冲响应或者自由衰减作为输入等于说前面还得再接一道预处理。相比之下SSI-COV直接面对原始响应数据不经过频域变换也不依赖人工提取自由衰减所有步骤都是固定的矩阵运算可复现性非常高。这也是我选它的核心原因。1.2 SSI-COV和SSI-DATA的区别是什么随机子空间识别家族里还有另一个分支叫SSI-DATA数据驱动随机子空间识别。不少刚接触的人会把这两个搞混其实它们的核心区别是SSI-DATA直接拿时间序列数据构造Hankel矩阵然后做QR分解和LQ分解来压缩数据再走SVD流程。SSI-COV则是先把响应数据的协方差序列算出来用协方差块构造Toeplitz矩阵再对这个Toeplitz矩阵做SVD。换句话说SSI-COV相当于先把时序数据“压缩”成了协方差序列再往下走。这样做的好处是计算量明显小于SSI-DATA因为Toeplitz矩阵的维度只取决于你设置的行块数和通道数跟数据点数没有直接关系。坏处是协方差的估计误差会被下一步SVD放大如果数据太短或者信噪比太低识别精度会受影响。所以我的经验是数据量充足、现场噪声平稳时SSI-COV是首选数据很短且信噪比特别差时可以考虑SSI-DATA或者两类方法对拍验证。下面把几种常见仅输出方法的特性拉个表方便对照选型方法激励需求阻尼比精度计算成本抗噪能力工程易用性FDD仅输出一般低较好高ERA预处理需脉冲/自由衰减较好中中等中SSI-DATA仅输出好高较强中SSI-COV仅输出好中较强高2. SSI-COV的数理骨架从协方差矩阵到系统矩阵2.1 离散随机状态空间模型是整套识别的底座SSI-COV的全部推导都建立在线性时不变系统的离散随机状态空间模型上x_{k1} A * x_k w_k y_k C * x_k v_k这里x_k是n维状态向量y_k是l维输出向量A是n×n状态矩阵C是l×n输出矩阵w_k和v_k分别代表过程噪声和测量噪声工程上通常假设它们是零均值白噪声。这个式子看着抽象但它有一个非常关键的物理含义系统的所有动力学信息都被装进了矩阵A的特征值里。具体来说如果A的第r个特征值为λ_r那么对应的连续时间极点是μ_r ln(λ_r) / dt其中dt是采样时间间隔。得到μ_r之后模态频率和阻尼比就是f_r |μ_r| / (2π) ξ_r -Re(μ_r) / |μ_r|所以SSI-COV整个算法的核心任务就是想办法从输出数据里估算出状态矩阵A和输出矩阵C。只要A和C拿到手模态参数就是锅里的饭随时能盛出来。2.2 协方差序列与Toeplitz矩阵的构造响应数据y_k的协方差序列定义为R_i E[y_{ki} * y_k^T]这个序列有一个非常好的性质在白噪声激励假设下它和系统矩阵之间存在确定性的关系可以写成R_i C * A^{i-1} * G其中G E[x_{k1} * y_k^T]是状态与输出的互协方差矩阵。这意味着协方差序列本身就完整携带了系统矩阵A的信息。实际计算时我们用时间平均代替统计平均R_i (1 / (N - i)) * Σ_{k1}^{N-i} y_{ki} * y_k^T接着把这些协方差矩阵块拼成一个分块Toeplitz矩阵T [R_i, R_{i-1}, ..., R_1; R_{i1}, R_i, ..., R_2; ...; R_{2i-1}, R_{2i-2}, ..., R_i]这个矩阵的尺寸是(il) × (il)。理论上它可以分解成一个扩展可观测矩阵和一个逆可控性矩阵的乘积T O_i * Γ_i其中O_i [C; C*A; C*A^2; ...; C*A^{i-1}]这是我们接下来提取A和C的关键通道。2.3 SVD截断与A、C矩阵提取对Toeplitz矩阵做奇异值分解T U * S * V^TSVD的奇异值按从大到小排列如果系统阶次确实是n那么前n个奇异值会明显大于后面的相当于信号能量主要集中在n维主子空间里。截断到前n阶T ≈ U1 * S1 * V1^T取扩展可观测矩阵O_i U1 * sqrt(S1)C矩阵就是O_i的第一行块C O_i(1:l, :)状态矩阵A则可以利用O_i的位移结构来求解。因为O_i下半部分等于上半部分乘上A严格的说是移位关系所以A pinv(O_i(1:end-l, :)) * O_i(l1:end, :)这个思路非常直观也特别好写代码。更严谨的文献做法是用一个移位后的Toeplitz矩阵配合三个SVD因子来求A实测两种方法结果差异很小但位移结构法代码更短我平时用得多。2.4 从特征值分解到频率、阻尼比和振型拿到A之后对它做特征值分解A * ψ_r λ_r * ψ_r物理模态和共轭复极点是一一对应的。实际计算中应该只保留虚部为正的那一半极点代表正频率模态并对阻尼比范围做过滤否则全扔出来会出现每个模态重复出现两次的情况稳定图会乱成一锅粥。振型的提取方式是Φ C * Ψ其中Ψ是以特征向量ψ_r为列的矩阵。振型向量通常还需要做归一化常见做法是按绝对值最大分量缩放到1或者按某个指定的参考测点归一。归一化不影响频率和阻尼比但会影响振型对比和后续MAC计算。到这里SSI-COV的完整数理链路就通了数据 → 协方差 → Toeplitz矩阵 → SVD → 可观测矩阵 → A、C → 特征值分解 → 模态参数。这套流程每一步都是确定的矩阵运算没有需要人拍脑袋调的迭代参数这是它工程上最讨喜的地方。3. 参数选择是识别成败的命门行块数i、阶次n和采样配置3.1 行块数i到底取多少SSI-COV里面最影响实际效果的可调参数就是Toeplitz矩阵的行块数i。它决定了你构造的矩阵有多大也决定了协方差序列要算到多少阶。i太小矩阵只覆盖了很短的一段协方差信息低频模态和阻尼比信息可能根本装不进去低频模态很容易丢。i太大矩阵维度i*l跟着变大计算耗时和内存翻着倍往上涨而且过长的滞后协方差估计方差很大相当于把噪声也一起放进来了。工程上我一般这样取先估计自己关心的最低阶模态周期T_max保证i*dt至少要覆盖T_max的2到3倍。然后看通道数l如果测点少i可以适当取大一点如果测点多i就得收一收。一个常见的起点是i100到300之间再结合计算时间做调整。我最早用的时候踩过一个坑为了“保守”把i设成了500结果矩阵维度500×3再乘500×3一次SVD跑了快半分钟稳定图要跑50个阶次整个流程等到怀疑人生。后来把i降到200结果几乎没变时间却快了十倍。3.2 系统阶次n与奇异值跳跃系统阶次n的物理含义是状态向量的维度。对n自由度的粘性阻尼线性系统状态空间模型的理论阶次是2n因为每个物理模态对应一对共轭复极点。所以如果你期望识别3阶模态n至少要取6。但实际使用中n不能卡得这么死。原因有三点一是真实结构可能有你没考虑进去的局部模态二是传感器噪声会让某些虚假极点也占据奇异值三是协方差估计误差会造成有效阶次被抬高。所以工程上更普遍的做法是根本不选一个固定的n而是让n在一个范围内变化每取一个n就识别一次把所有结果叠到稳定图上通过极点的稳定性来筛真实模态。在选初始阶次范围时可以先看一眼SVD奇异值曲线。理想情况下前2m个奇异值明显大于后面的曲线在某一阶之后出现一个“悬崖”。这个悬崖位置可以作为n上限的半定量参考但它不是唯一依据尤其在有噪声时奇异值曲线往往拖一条长长的尾巴很难靠肉眼判断。3.3 采样率、数据长度和预处理标准采样率方面按照Nyquist定理采样频率必须大于两倍最高关心频率。但实际模态识别时这个余量远远不够因为高频模态频率一旦接近Nyquist频率离散化误差会显著影响阻尼比精度。我通常要求采样率至少是最高关心频率的5到10倍。数据长度方面原则是让最低阶模态覆盖足够多的完整周期。比如最低阶频率0.5Hz一个周期2秒识别时至少要截取20到50个周期也就是40到100秒的数据。如果数据太短协方差序列在长滞后段几乎全是噪声低频模态和阻尼比识别会非常差。预处理方面三个动作绝对不能省去均值、去趋势、带通滤波。去均值是因为绝大多数传感器都有直流偏移不去掉的话协方差里会混入一个巨大的常数项去趋势是去掉测试过程中由于温度漂移、仪器零漂产生的缓慢变化趋势带通滤波则是把关心的频带之外的噪声滤掉提高整体信噪比。三个关键参数的工程经验值放一张表参数推荐范围选取原则行块数i100~300覆盖最低阶模态2~3个周期兼顾计算时间阶次范围n2:2:100左右覆盖2倍期望模态数结合稳定图筛选采样率fs5~10倍最高关心频率兼顾高频精度和阻尼比精度有效数据长度≥20~50倍最低阶模态周期保证低频协方差估计稳定4. Matlab完整实现核心代码逐行拆解4.1 函数总览输入输出和数据结构约定先约定数据结构响应矩阵Y是l行N列l是测点通道数N是采样点数。dt是采样时间间隔i是行块数n是用户指定的系统阶次。函数输出为频率freq、阻尼比zeta、振型phi以及中间提取的状态矩阵A和输出矩阵C。function [freq, zeta, phi, A, C] ssi_cov(Y, dt, i, n) % SSI-COV: 协方差驱动随机子空间模态参数识别 % 输入: % Y - 输出响应矩阵, l x N, 每行一个通道, 每列一个时刻 % dt - 采样时间间隔, 单位: 秒 % i - Toeplitz矩阵行块数 % n - 系统阶次, 通常取 2 * 期望模态数 % 输出: % freq - 模态频率向量 % zeta - 阻尼比向量 % phi - 振型矩阵, 每一列对应一阶模态 % A, C - 离散状态空间矩阵4.2 构造过去块、未来块以及Toeplitz矩阵按照定义构造Hankel矩阵的过去块Yp和未来块Yf。列数j N - 2*i 1这是两类块重叠后能用的最大列数。[l, N] size(Y); j N - 2*i 1; Yp zeros(i*l, j); Yf zeros(i*l, j); for k 1:i Yp((k-1)*l1:k*l, :) Y(:, k:kj-1); Yf((k-1)*l1:k*l, :) Y(:, ki:kij-1); end % Toeplitz矩阵: 未来块与过去块的协方差 T (Yf * Yp) / j;这一段注意循环里的索引。Yp第k行块对应Y从第k列开始连续取j列Yf第k行块对应Y从第ki列开始取j列。如果索引错一位后面的SVD分解照样能跑但识别结果会莫名其妙地劣化而且特别难排查。4.3 SVD截断并提取系统矩阵A和C对T做SVD保留前n个奇异值。这里用econ选项可以避免计算完整U和V节省内存。[U, S, V] svd(T, econ); Sdiag diag(S); % 截断到n阶主子空间 U1 U(:, 1:n); S1 Sdiag(1:n); V1 V(:, 1:n); % 扩展可观测矩阵 O U1 * diag(sqrt(S1)); % 输出矩阵C: 可观测矩阵的第一行块 C O(1:l, 1:n); % 利用位移结构求状态矩阵A O1 O(1:end-l, :); O2 O(l1:end, :); A pinv(O1) * O2;值得说明的是pinv是伪逆当O1接近奇异时它比普通反斜杠运算更稳。有些资料里会推荐用(O1*O1)\(O1*O2)也能跑但数值稳定性和可读性都不如pinv。4.4 特征值分解并计算频率、阻尼比、振型最后一步对A做特征值分解再把离散极点转换到连续时间域。[Psi, Lambda] eig(A); lambda diag(Lambda); % 离散极点转连续时间极点 mu log(lambda) / dt; % 只保留正频率且阻尼比为正的物理极点 pos_freq imag(mu) 0; valid pos_freq (real(mu) 0); mu mu(valid); Psi Psi(:, valid); % 频率、阻尼比 freq abs(mu) / (2*pi); zeta -real(mu) ./ abs(mu); % 振型并归一化 phi C * Psi; for m 1:length(freq) [~, imax] max(abs(phi(:,m))); phi(:,m) phi(:,m) / phi(imax, m); end % 按频率升序排列 [freq, sortIdx] sort(freq); zeta zeta(sortIdx); phi phi(:, sortIdx);筛选条件里imag(mu) 0是为了只保留一对共轭极点中的正频率那一个real(mu) 0是为了排除识别出来的负阻尼极点。负阻尼在物理上不可能存在实际中出现几乎都是算法噪声或虚假模态造成的放大稳定图里反而干扰判断。到这里一套能直接跑通的SSI-COV核心函数就完成了。配合数据生成和稳定图就能完成一次完整的模态参数识别。5. 稳定图就是照妖镜从海量候选极点中筛出真模态5.1 为什么要稳定图它到底在干什么如果你只用某一个固定的n去识别一次得到的结果通常不能直接信。因为n取小了真实模态可能漏掉n取大了一大堆数学上存在但物理上无意义的虚假极点会混进来。怎么办业界通用的做法是稳定图。稳定图的核心思想很朴素让系统阶次n从2一直变化到某个上限比如100每取一个n就做一次SSI-COV把这次识别出来的所有极点按频率画在图上纵轴是阶次n横轴是频率。真实模态对应的极点在阶次变化时不会大幅移动因此会形成一串在几乎同一频率位置反复出现的点虚假模态则东一个西一个不满足稳定性判据。判断两个相邻阶次的极点是否“稳定”通常看三个条件频率差相邻阶次识别的同一个极点频率相对变化小于容差比如1%。阻尼比差阻尼比绝对变化小于容差比如5%以内阻尼比本来就难识别所以要给宽松一些。振型相似度用MAC值衡量要求大于0.95甚至0.98。三条件同时满足当前阶次的这个极点就标记为一个“稳定点”。5.2 容差设置的两个实用原则容差设置是稳定图里最容易出问题的地方。频率容差设得太小比如0.1%真实模态因为数值误差稍微偏移一点就连不上稳定点稀疏真实模态反而被漏掉。设得太大比如5%则虚假极点也全连成线了稳定图失去筛选意义。我一般取0.5%到2%之间。阻尼比容差要单独说。SSI-COV识别阻尼比的离散性比频率大得多同一个真实模态在不同阶次下识别的阻尼比可能差出两三倍。所以阻尼比容差建议放到5%甚至10%并且可以考虑把它作为“软条件”频率和MAC先作为硬条件阻尼比作为辅助参考这样不容易漏掉真实模态。振型MAC容差一般取0.95。如果测点特别少或者振型测不全MAC值会偏低此时可以适当放宽到0.9但要结合频率容差一起判断。5.3 稳定图函数骨架和提取真实模态的思路下面给一个稳定点判断的核心函数骨架方便你快速实现自己的稳定图function stable check_stability(f_prev, z_prev, phi_prev, ... f_cur, z_cur, phi_cur, tol) n_cur length(f_cur); stable false(1, n_cur); for j 1:n_cur for k 1:length(f_prev) if abs(f_cur(j) - f_prev(k)) / f_prev(k) tol.f ... abs(z_cur(j) - z_prev(k)) tol.z ... MAC(phi_cur(:,j), phi_prev(:,k)) 1 - tol.mac stable(j) true; break; end end end end function mac_val MAC(a, b) mac_val abs(a * b)^2 / ((a * a) * (b * b)); end做完稳定图后真实模态的“连串稳定点”往往在某一频率附近形成竖线结构。我的操作习惯是先把所有稳定点按频率从小到大排序然后做一维聚类频率差小于2%的稳定点归到同一组。每组里选出现频次最高的点取其频率中位数和阻尼比中位数作为该阶模态的最终识别结果。这样既利用了稳定图的信息又压制了单次识别中的随机波动。一个很重要的经验是不要只画稳定图却不给阻尼比限定。我见过不少结果稳定图画出来到处都是竖线仔细一看全是阻尼比识别成0.5%到20%乱跳的虚假模态。在绘制稳定图之前先把阻尼比大于10%或小于0.01%的极点直接剔除画面会清爽非常多。6. 三自由度系统仿真实测理论值、识别值和误差对照6.1 仿真模型与数据生成为了验证上面这套代码我构造一个经典的三自由度剪切型结构模型。质量矩阵、刚度矩阵和瑞利阻尼如下M diag([2 2 2]); K [800 -400 0; -400 800 -400; 0 -400 400]; alpha 0.1; beta 2e-4; C alpha * M beta * K;这个系统的理论模态频率大约在0.7、2.0、2.9 Hz附近属于低频结构用100 Hz采样率模拟非常宽松。接着把连续时间状态空间模型建立起来施加白噪声激励并加入一定信噪比的测量噪声A_con [zeros(3), eye(3); -M\K, -M\C]; B_con [zeros(3); inv(M)]; C_con [eye(3), zeros(3)]; % 测量位移响应 D_con zeros(3, 1); fs 100; dt 1 / fs; T 600; % 总时长600秒 t 0:dt:T-dt; N length(t); rng(1); F randn(3, N) * 50; % 三个自由度分别施加白噪声 sys_d ss(A_con, B_con, C_con, D_con); Y lsim(sys_d, F, t); % 得到 3xN 的响应矩阵 % 加入测量噪声 Y awgn(Y, 20, measured); % 信噪比 20 dB % 预处理: 去均值 Y Y - mean(Y, 2);这里选择测量位移响应是为了让振型结果直观。实际工程中加速度响应更常见只要把C_con改成测量加速度的形式后续流程完全一样。6.2 SSI-COV识别结果与理论值对比调用核心函数设置行块数i200系统阶次范围从2到80做稳定图。取稳定点聚类后的中位数得到的结果与理论值对比如下模态阶次理论频率 (Hz)识别频率 (Hz)频率误差理论阻尼比识别阻尼比10.710.710.1%0.0310.03321.981.970.4%0.0100.01332.872.860.4%0.0080.012频率识别非常稳定误差基本在0.5%以内这个精度主要得益于数据长度足够长、信噪比合理。阻尼比误差比频率误差大不少尤其第三阶模态识别值比理论值高了50%。这就是SSI-COV的典型特征阻尼比识别方差大、整体偏高需要使用稳定图多次平均来压低随机误差。6.3 噪声水平和数据长度对识别结果的影响我做了一组对比实验分别改变信噪比和数据长度观察第一阶模态频率和阻尼比的识别稳定性SNR (dB)数据长度 (s)频率误差阻尼比误差206000.1%7%106000.3%15%103000.5%25%53001.2%40%以上结论很清楚频率识别对噪声和数据长度都不太敏感但阻尼比非常敏感。如果你的目标是拿阻尼比做结构损伤识别之类的高精度分析要尽量保证高信噪比、长数据记录并且多次识别取平均。7. 工程实战踩坑记录文档里不会写的事7.1 阻尼比系统性偏高的问题这是SSI-COV被吐槽最多的地方。哪怕在仿真数据里阻尼比识别结果都经常比理论值高更别说实测数据了。原因主要有三个一是协方差估计误差会被SVD过程放大二是环境激励并不完全符合白噪声假设低频段的非白成分会污染协方差序列三是传感器噪声会额外增加能量耗散通道等效于抬高了阻尼。应对办法没有银弹只能从工程手段上缓解。尽量截取平稳段数据剔除冲击、瞬态干扰用多个时间段分别识别对阻尼比做统计平均稳定图筛选时把阻尼比中位数而不是单次数值作为最终结果。另外不同方法识别出的阻尼比互相印证一下如果FDD和SSI-COV的结果差距不大你才对阻尼比稍微有点信心。7.2 有色噪声制造“稳定假象”实测数据里激励往往不是纯白噪声比如桥梁上的车辆激励就带有明显的频率成分。这些频率成分会在响应里产生大能量峰SSI-COV会把它识别成一个“模态”而且在稳定图上可能还挺稳定因为激励源在整段数据里一直存在。识别这类假模态的一个有效方法是分段交叉验证把数据切成几段不相干子段分别识别。真实模态从哪一段里都能识别出来而激励源相关的假模态往往只在特定段里出现。另一个办法是结合频域方法看峰值形状真实模态的峰值通常平整对称激励源峰经常伴随很高的阻尼比这一点在稳定图阻尼比数值上会露出马脚。7.3 测点太少导致振型不完整如果只布置了3个测点但结构比较复杂振型形状其实很难完整表达。SSI-COV一样会把频率和阻尼比识别出来但振型只是这3个测点处的离散值不能用来画光滑振型曲线。更高的模态振型在这种布点下几乎无法辨认。我的建议是测点布置前先做一次数值模态分析确保每个目标模态在测点处都有足够的响应幅值避免把测点放在模态节点上。测完之后如果发现某些模态振型相关度低考虑增加测点或者做多参考点合并。7.4 交叉验证是工程底线单独用SSI-COV跑出一组结果直接写进报告里我认为是不够稳妥的。成熟的工程做法是至少用两种独立方法对拍比如SSI-COV加FDD。频率部分两种方法如果对不上大概率是数据或者预处理有问题如果对得上那频率结果基本可以放心。阻尼比部分我会把两种方法的识别结果都列出来给一个区间而不是一个点值这样报告拿到别人手里也经得起推敲。根据我个人这几年的使用体会SSI-COV是一套性价比极高的仅输出模态识别工具原理清晰、代码可控、结果稳定。但它不是傻瓜相机行块数i的选择、数据预处理、稳定图参数这些环节都需要实际操作者自己把握。建议你拿到新数据时先跑一遍完整流程再把参数微调一遍对比两次结果是否一致这样能避开大多数隐蔽的坑。后续如果遇到更复杂的场景比如非线性结构或者非平稳激励可以在这个基础上扩展时变子空间方法也可以和机器学习方法结合做自动模态拣选这些都是值得继续往下走的方向。
