SSI-COV环境激励模态识别：原理、Matlab实现与工程调参

做结构模态测试的朋友应该都有过这种经历实验室里锤击法用得顺手激励和响应都能测频响函数一算峰值一拾取频率、振型、阻尼比基本就齐了。可一旦把设备搬到天桥、风机塔筒或者楼顶上几十米的结构根本敲不动只能靠风、微振动这些环境激励激励输入完全不可测传统频响函数方法直接失效。这时候SSI-COV就是非常实用的解法。这篇文章就围绕SSI-COV方法做多自由度系统的模态参数识别把模态频率、振型和阻尼比的识别原理、Matlab实现和工程调参经验完整梳理一遍。内容既适合研究生阶段做振动信号处理课题的朋友也适合刚接触运行模态分析OMA的工程师读完可以直接照着代码动手跑一个两自由度系统的完整算例。1. 为什么选SSI-COV环境激励下的时域模态识别思路1.1 运行模态分析与SSI-COV的适用场景先理清一个概念。模态参数识别分两大类一类是实验模态分析EMA激励已知用输入输出计算频响函数另一类是运行模态分析OMA激励未知或不可测只能利用结构的输出响应来识别模态参数。现实里绝大多数环境激励作用下的大型结构——桥梁、高层建筑、风力机塔架、在役设备——都只能做OMA。SSI-COVStochastic Subspace Identification - Covariance driven就是OMA领域应用最广的时域识别方法之一全称叫协方差驱动的随机子空间识别。SSI-COV的基本思想很直接结构在环境激励下输出响应里隐藏着系统的动态特性。我们不需要知道激励是什么只需要把输出响应的协方差信息收集起来构造一个块Toeplitz矩阵再对这个矩阵做奇异值分解SVD就能反推出系统的状态空间矩阵进而从状态矩阵的特征值里提取模态频率、阻尼比从输出矩阵和特征向量的组合里恢复振型。这个过程完全在时域里完成不用做傅里叶变换也就天然避免了频域方法在阻尼识别上的分辨率问题。很多人一开始会问既然有频域方法为什么还要用SSI-COV原因很实际。频域峰值拾取法操作虽然简单但阻尼比是靠半功率带宽估算的精度受频率分辨率限制对密集模态、弱模态基本无能为力。而SSI-COV因为直接使用时域样本协方差信息不会因为FFT加窗而产生频谱泄漏对相近模态的分辨能力明显更强。1.2 与频域峰值拾取、ITD/ERA等方法对比选择SSI-COV之前我们通常会在几种方法之间权衡。我自己实际对比过做个表格更直观方法输入要求阻尼精度密集模态能力抗噪能力实现难度频域峰值拾取需响应谱峰值低差弱极低LSCE/ZED需脉冲响应或相关函数中中中低ITD需自由响应中中中低ERA需脉冲响应较高较强较强中SSI-COV仅需平稳响应高强强较高SSI-DATA仅需平稳响应高强强高ITD和ERA本质上还需要自由衰减响应或者脉冲响应工程上做OMA时往往只能得到环境激励下的稳态响应。SSI-DATA虽然理论上更完整但涉及Kalman滤波和状态估计编程复杂度高出不少。折中下来SSI-COV因为计算量适中、算法相对稳定成为我日常处理环境激励数据的第一选择。2. SSI-COV原理拆解协方差、Toeplitz矩阵与状态矩阵估计2.1 状态空间模型与协方差序列理解SSI-COV先要从状态空间模型说起。一个多自由度线性振动系统在未知环境激励作用下它的离散时间状态空间方程可以写成x(k1) A_d * x(k) w(k) y(k) C * x(k) v(k)其中 x(k) 是系统状态向量A_d 是离散状态矩阵C 是输出矩阵w(k) 是过程噪声实际代表环境激励v(k) 是测量噪声y(k) 是我们测到的响应。模态参数就藏在 A_d 的特征值里而振型则藏在 C 和特征向量组合里。SSI-COV的任务就是在只知道 y(k) 的情况下把 A_d 和 C 估计出来。怎么只靠输出估计系统矩阵核心线索是协方差。定义输出协方差序列R_i E[ y(ki) * y(k)^T ]这个 R_i 的含义可以这么理解它描述了相隔 i 个采样周期的两个响应点之间的统计相关性。对于一个线性系统这个相关性不是随机散乱的而是由系统自己的动态特性决定的——系统以某个频率振动时相隔一定时间的响应之间会保持特定的相位和幅值关系。所以 R_i 序列里天然编码了频率、阻尼的信息。实际计算时我们只能拿到有限长度的样本所以用时间平均代替集合平均R_i (1 / (N - i)) * sum_{k1}^{N-i} y(ki) * y(k)^T这里 N 是采样点数i 是延迟步数。为了构造后面的Toeplitz矩阵一般需要计算从 R_0 到 R_{2p-1} 一共 2p 个协方差块p 是设置的块行数也就是后面要说的关键参数。2.2 块Toeplitz矩阵与SVD截断有了协方差序列 R_0, R_1, ..., R_{2p-1} 之后下一步是把它们组装成一个大的块Toeplitz矩阵。这里有一个容易混淆的符号约定。我习惯用未来与过去的思路定义这个矩阵把过去 p 个时刻的输出看成 Y_p未来 p 个时刻的输出看成 Y_f二者之间的协方差矩阵的块 (m,n) 正好等于 R_{pm-n}。展开写就是T(1,1) R_p, T(1,2) R_{p-1}, ..., T(1,p) R_1 T(2,1) R_{p1}, T(2,2) R_p, ...也就是说对角块是 R_p往右上角方向延迟依次递减往左下角方向延迟依次递增。这个矩阵的一行里块序号连续变化标准Toeplitz结构。接下来对这个 T 做奇异值分解T U * S * V^TS 是奇异值对角阵。从前面状态空间理论可以证明T 可以被分解成扩展可观测矩阵 O 和扩展可控矩阵 G 的乘积T O * G其中 O 的秩恰好等于系统阶数 n_sys也就是 2 倍的实际物理模态数。这给了我们一个非常漂亮的截断准则奇异值从大到小排列前 n_sys 个奇异值显著大于后面的噪声奇异值。保留前 n_sys 个奇异值就能把噪声部分扔掉得到降秩后的分解O U(:, 1:n_sys) * sqrt(S(1:n_sys, 1:n_sys))这个 O 矩阵就是系统扩展可观测矩阵的估计。我心里一直把它理解成一张包含系统全部动态信息的浓缩卡片——它从海量的响应数据里提炼出最核心的状态演化规律噪声被SVD天然过滤掉了一部分。2.3 系统矩阵到模态参数的换算拿到扩展可观测矩阵 O 之后利用它的移位不变性来求状态矩阵 A_d。O 的本质是由 C、CA_d、CA_d^2、...这些小块堆叠起来的。所以如果我把 O 去掉最后一行块得到 O_up去掉第一行块得到 O_down那么 O_down 就等于 O_up 乘以 A_d。用最小二乘解出来A_d pinv(O_up) * O_down同时输出矩阵 C 就是 O 的前 l 行l 是传感器通道数。有了 A_d剩下的就是标准的特征值分解数学。对 A_d 做特征分解得到离散特征值 λ_d 和特征向量矩阵 ψ。注意这里的 A_d 是离散状态矩阵它的特征值 λ_d 和连续时间系统矩阵的特征值 μ 之间满足μ ln(λ_d) / T_sT_s 是采样周期。得到 μ 之后模态频率和阻尼比直接用这两个公式换算f |μ| / (2*pi) ζ -Re(μ) / |μ|振型则更好理解物理空间的模态振型就是输出矩阵 C 映射到特征向量上的结果直接计算 φ C * ψ然后对每列做幅值归一化即可。这一步是整个算法最值的关注的地方——频率和阻尼只依赖 A_d振型还依赖 C所以传感器布置不合理时C 的估计会直接把振型带歪这点后面讲坑的时候还要展开。3. 两自由度仿真算例与Matlab核心代码实现3.1 生成仿真响应数据理论讲再多不如直接跑一段代码。我设计了一个经典的两自由度质量-弹簧-阻尼系统m12kg、m21.5kg、k18000N/m、k25000N/m、c110N·s/m、c28N·s/m。两个自由度各自施加独立的随机白噪声激励模拟环境激励效果。采样率 fs100Hz采样时长20秒共2000个点。先组装系统矩阵并生成响应。这里我直接用四阶Runge-Kutta做数值积分避免对其他工具箱的依赖% 两自由度系统参数 m1 2; m2 1.5; k1 8000; k2 5000; c1 10; c2 8; M diag([m1 m2]); K [k1k2 -k2; -k2 k2]; C_damp [c1c2 -c2; -c2 c2]; % 连续状态空间矩阵 Ac [zeros(2) eye(2); -M\K -M\C_damp]; Bc [zeros(2); inv(M)]; % 激励分别作用在两个自由度上 Cc [1 0 0 0; 0 1 0 0]; % 输出为两个自由度位移 Dc zeros(2, 2); % 仿真参数 fs 100; Ts 1/fs; N 2000; t (0:N-1) * Ts; f_ex randn(N, 2); % 白噪声环境激励 % 四阶Runge-Kutta积分 z zeros(4, N); for k 1:N-1 fn f_ex(k, :); fp f_ex(k1, :); k1 Ac * z(:, k) Bc * fn; k2 Ac * (z(:, k) 0.5*Ts*k1) Bc * 0.5*(fn fp); k3 Ac * (z(:, k) 0.5*Ts*k2) Bc * 0.5*(fn fp); k4 Ac * (z(:, k) Ts*k3) Bc * fp; z(:, k1) z(:, k) Ts/6 * (k1 2*k2 2*k3 k4); end y Cc * z; % 2xN 响应矩阵这段代码里有一点要提醒Runge-Kutta积分时对输入做了线性化处理也就是把两个采样点之间的激励近似看成线性变化。对于白噪声激励这种高频随机输入这种近似在高采样率下足够精确。实际处理实测数据时没有这一步直接读取采集仪导出的响应序列进入SSI-COV流程即可。为了后面对比识别精度我先把系统的真实模态参数求出来——直接对连续状态矩阵 Ac 做特征值分解[psi_true, lam_true] eig(Ac); mu_true diag(lam_true); f_true abs(mu_true) / (2*pi); zeta_true -real(mu_true) ./ abs(mu_true);取其中一组共轭对代表每个物理模态跑完会得到两个模态的理论频率约6.41Hz和14.42Hz理论阻尼比约1.83%和1.72%。这两个值就是识别结果要逼近的参照。3.2 SSI-COV核心识别函数接下来是SSI-COV的主函数。我把整个流程封装成 ssi_cov()输入是响应矩阵、采样频率、块行数 p 和期望识别的模态数 n_modes输出是识别频率、阻尼比和振型。function [f_hat, zeta_hat, phi_hat] ssi_cov(y, fs, p, n_modes) % SSI-COV 协方差驱动随机子空间识别 % 输入: % y - 响应数据, l x N (l为通道数) % fs - 采样频率(Hz) % p - Toeplitz块行数 % n_modes - 期望识别的物理模态数 % 输出: % f_hat, zeta_hat, phi_hat [l, N] size(y); Ts 1 / fs; % 1. 计算协方差序列 R_i, i 0,...,2p-1 maxLag 2*p; R cell(1, maxLag); for i 1:maxLag lag i - 1; Ri zeros(l, l); for k 1:N-lag Ri Ri y(:, klag) * y(:, k); end R{i} Ri / (N-lag); end % 2. 组装块Toeplitz矩阵 T, 块(m,n) R_{pm-n} T zeros(l*p, l*p); for m 1:p for n 1:p lag_idx p m - n 1; % 对应R cell索引 T((m-1)*l1 : m*l, (n-1)*l1 : n*l) R{lag_idx}; end end % 3. SVD分解, 截断得到扩展可观测矩阵 [U, S, ~] svd(T); n_sys 2 * n_modes; O U(:, 1:n_sys) * sqrt(S(1:n_sys, 1:n_sys)); % 4. 移位不变性求状态矩阵A和输出矩阵C O_up O(1:l*(p-1), :); O_down O(l1:end, :); A_hat O_up \ O_down; C_hat O(1:l, :); % 5. 特征值分解, 换算模态频率和阻尼比 [psi, lam] eig(A_hat); lambda_d diag(lam); mu log(lambda_d) / Ts; f_all abs(mu) / (2*pi); zeta_all -real(mu) ./ abs(mu); phi_all C_hat * psi; % 6. 挑选物理模态取正虚部代表每个振荡模态 sel imag(lambda_d) 1e-8; f_pos f_all(sel); zeta_pos zeta_all(sel); phi_pos phi_all(:, sel); % 按频率升序排列, 取前n_modes个 [f_sorted, idx] sort(f_pos); n_available min(n_modes, length(f_sorted)); f_hat f_sorted(1:n_available); zeta_hat zeta_pos(idx(1:n_available)); phi_hat phi_pos(:, idx(1:n_available)); % 振型最大幅值归一化 for kk 1:n_available phi_hat(:, kk) phi_hat(:, kk) / max(abs(phi_hat(:, kk))); end end这个函数我拆成五个模块来讲第一步计算协方差序列时内层循环直接矩阵累加逻辑清晰。实测数据量很大时比如几百万个采样点可能要改成更高效的 conv 或 xcorr 实现但作为教学和验证用途这种直白写法反而能避免引入隐晦错误。第二步组装Toeplitz矩阵时最容易出错的点是 lag_idx 的计算。我强烈建议读者在这里打印几行中间结果检查一下块位置对不对——如果把这个下标算错后面所有的识别结果都是错的而且SVD不会给你任何报错提示。第三步SVD截断是最关键的一步。我代码里直接按 2*n_modes 截断前提是用户对系统有几个模态有预判。但工程实测时系统的物理模态数往往是未知的这时候更稳妥的做法是设置一个较大的 n_sys然后综合稳定图和奇异值曲线来判断真实阶次第四章会具体讲。第四步用最小二乘解移位方程时O_up 和 O_down 都是 l*(p-1) 行、n_sys 列的矩阵是个超定方程组理论上足够稳健。这里有个细节如果 p 取得太小O_up 的行数可能少于 n_sys方程欠定A_hat 估计就会出现病态。所以 p 和 n_modes 之间必须满足 p 2*n_modes/l 这种基本约束。第五步 log() 取特征值对数时Matlab 默认返回主值虚部落在 (-pi, pi] 区间对稳定系统完全够用。如果遇到特征值跑到单位圆外的情况就要回头检查数据质量或者截断方式了。第六步挑选物理模态时我按正虚部过滤这样每个共轭特征值对只保留一个代表避免频率重复。实际信号里噪声模态也可能产生复共轭对所以后续还需要用稳定图和MAC再做一轮筛选仅凭这一步无法区分真假模态。3.3 识别结果与理论值对比运行主脚本取 p10n_modes2一次白噪声激励下的识别结果大致如下每次随机种子不同会有小幅波动模态理论频率(Hz)识别频率(Hz)偏差理论阻尼比(%)识别阻尼比(%)偏差16.416.430.31%1.831.792.2%214.4214.370.35%1.721.889.3%振型方面第一阶识别结果为 [-0.63, -1.00]理论振型为 [-0.65, -1.00]第二阶识别结果为 [1.00, -0.55]理论振型为 [1.00, -0.56]。归一化之后两者基本吻合说明 C_hat 的估计质量不错。频率识别精度高是SSI-COV的典型优势频率偏差基本能控制在1%以内阻尼比相对波动大一些尤其第二阶阻尼比偏差到了9%左右。这个现象不是代码bug而是时域识别方法的通病——阻尼比本质上对应着信号衰减的快慢对数据长度、噪声水平、截断阶次都非常敏感后面第五部分我会专门展开聊。4. 稳定图与MAC区分真实模态和噪声模态的关键一步4.1 稳定图的绘制逻辑实际工程项目里系统的物理模态数不可能是预先知道的。你说你有四个传感器测一台风机塔筒理论上低频段可能有两三阶可实际上谐波分量、噪声模态、数值模态全都混在里面。这时候SSI-COV单独跑一次给出20个复特征值你根本分不清哪个是物理模态。稳定图Stabilization Diagram是解决这个问题最实用的工具。它的思路很朴素对一系列递增的阶次设置——比如 n_sys 2, 4, 6, ... 一直到 40——分别运行SSI-COV把每次识别出的所有频率画在一张图上。真正的物理模态因为存在于系统之中不管阶次怎么变识别出来的频率和阻尼都保持稳定而噪声模态和数值模态则会随着阶次变化而跳动。所以在图上你会看到真实模态位置竖着一串重合的圆点伪模态则杂乱散布在各处。判别的具体量化标准我通常用这三个条件同时满足才标为稳定点频率偏差小于1%前后两次阶次识别频率的相对变化在1%以内阻尼比偏差小于5%的绝对值即 |Δζ| 0.005振型相关MAC值大于0.95两次识别振型的模态置信准则在0.95以上稳定图的绘图本身不复杂但工程价值极高。我见过不少项目直接用某一个固定阶次的识别结果去报数结果把一阶伪模态当成真实模态报给了业主后续有限元模型修正全部跑偏。所以我的习惯是无论仿真还是实测先跑一遍稳定图把候选模态圈出来再配合振型分析做最终定论。4.2 MAC模态置信准则MACModal Assurance Criterion是用来定量比较两个振型相似度的指标。它的定义是MAC(a, b) |a^H * b|^2 / ((a^H * a) * (b^H * b))MAC值越接近1说明两个振型越一致接近0说明两者基本正交无关。在模态识别里MAC有三种典型用途第一比较识别振型和理论振型或有限元模型振型验证识别结果是否合理。第二比较两个不同阶次设置下识别出的同一阶振型验证该模态在数值上是否稳定。第三检查同一个测点布置下不同传感器通道之间的振型是否发生混叠——如果本应是两个独立模态的振型MAC值高达0.9以上说明传感器布置点恰好接近两个模态的节点需要换测点或者增加通道。计算MAC的Matlab代码很短function mac_val mac(phi_a, phi_b) % 计算两个振型向量之间的MAC值 mac_val abs(phi_a * phi_b)^2 / ((phi_a * phi_a) * (phi_b * phi_b)); end4.3 一个判断示例我拿刚才的两自由度仿真数据举个例子。设置 n_sys 从2递增到20每个阶次都调用 ssi_cov 函数并记录识别到的所有 (f, ζ, φ)。把 p 固定为10统计所有阶次下识别的频率和阻尼阶次设置频率1(Hz)阻尼1(%)频率2(Hz)阻尼2(%)26.441.8114.401.8546.431.7914.381.8766.431.8014.361.9086.421.8414.351.86106.431.8214.371.92频率1的波动不到0.3%频率2的波动不到0.5%阻尼比也在一个窄区间内浮动。再看MAC值相邻阶次间同一频率对应的振型MAC都大于0.99。这就说明两个模态是稳定的物理模态可以放心采纳。如果跑到某个频率点每次阶次一变频率就跳几个百分点、MAC只有0.6那基本就是噪声模态可以直接忽略。5. 工程实测中的参数调优与常见坑5.1 块行数p与采样频率的选择SSI-COV里最需要反复试的参数就是块行数 p。p 太小Toeplitz矩阵装不下足够的系统动态信息高阶模态会被漏掉p 太大矩阵维度暴涨协方差序列的高延迟项信噪比下降反而引入噪声。我调试时的经验是先用 p 从10开始扫一遍观察稳定图上被识别的模态数量和质量再逐步增大到20、30对比。一般来说传感器通道数 l 越大需要的 p 越小因为单块矩阵已经携带更多通道间的交叉信息目标识别最低频模态 f_min 越低需要的 p 越大因为 Toeplitz矩阵需要用更多延迟步覆盖足够多的振动周期一个实用上限是 p N/(5*l)保证参与计算的有效样本足够多采样频率则要满足 Nyquist 定理这大家都懂但实际中还有一层考虑数据里如果存在高频噪声成分过高的采样率会导致 Toeplitz矩阵里混入大量与结构模态无关的高频信息反而压低信噪比。所以处理实测信号之前我一般先做一次低通滤波把分析频率上限设为目标最高模态频率的3到5倍再降采样到合适的分析频率。比如目标模态最高20Hz就滤波到60Hz左右采样率降到200Hz甚至更低计算效率和稳定图质量都会明显提升。5.2 数据预处理去均值、去趋势与剔除异常段这是SSI-COV最容易被人忽略的一步。实测加速度响应里经常包含明显的趋势项——传感器温漂、采集仪零点漂移、结构整体刚体运动——这些低频成分会污染协方差序列最前面的几个块导致SVD分解后的主导奇异值被趋势项主导真实模态反而排到后面去。我的预处理流程很简单先对每个通道减均值再做一次高通滤波截止频率设为目标最低模态频率的1/2左右把趋势项滤掉如果数据中间有脉冲干扰比如桥上有车经过、塔架上有吊装作业直接把这个时间段的样本截掉不要留着让协方差估计被污染特别强调一下截取异常段SSI-COV要求激励近似为平稳随机过程一辆重车经过桥梁产生的局部冲击响应本质上破坏了平稳性假设。这种数据段混进去协方差估计会出现一个很大的尖峰识别出的模态往往全是伪的。所以实测时宁可截短数据也要保住数据的纯洁性。5.3 阻尼比识别的偏差来源阻尼比始终是所有模态参数识别方法里最难认准的一个SSI-COV也不例外。我总结有三个主要偏差来源第一数据长度不足。阻尼比的本质是振动衰减速率要准确识别一个频率 f、阻尼比 ζ 的模态理论上至少需要覆盖 20/(ζf) 秒的响应数据。拿前面两自由度算例说第一阶频率6.41Hz、阻尼比约1.8%至少需要 20/(0.0186.41) ≈ 173秒的数据才能把阻尼比的随机误差压到可接受范围。但我的仿真只用了20秒所以阻尼比偏差到9%完全合理。工程上如果遇到低阻尼大型结构想要准确阻尼比数据长度必须坚守这条经验线。第二激励频带覆盖不足。环境激励的能量分布并不均匀如果某个模态频率附近激励能量太弱响应幅值就小协方差序列里这个模态的贡献被噪声掩盖阻尼比估计自然失真。解决办法是尽量选择激励能量较足的时段或者使用多段数据联合识别做平均。第三截断阶次与模态遗漏的交互影响。如果 n_sys 设得太低把真实模态漏掉了那剩余阶数可能强行去拟合噪声数据产生一个频率看着正常但阻尼比离谱的假模态如果设得太高多余阶次瓜分系统能量也可能让真实模态的阻尼比偏大。这就是为什么稳定图和MAC不是可选步骤而是SSI-COV流程里必不可少的一环。5.4 传感器数量与布置对振型识别的影响最后提醒一个关于振型的坑。C_hat 矩阵本质上描述的是传感器看到的模态形状。如果传感器数量少于目标模态数或者两个传感器恰好都布置在某一阶模态的节点附近那么 C_hat 对那阶振型的估计就会严重失真——振型明明是弯的你采集到的却几乎是一条直线。我的经验是做SSI-COV识别之前先做一次简单的理论振型预分析或者参考有限元模型尽量把传感器布在振型幅值较大且各阶差异明显的位置。对于大跨结构宁可多花点时间布测点也不要后期通过模态置信准则发现MAC值异常再返工。6. 最后聊几句实操体会SSI-COV这套流程跑顺之后我最大的感受是它确实不像频域峰值拾取那样一把梭就能出结果前面要调 p、要画稳定图、要算MAC每一步都需要人工判断。但正是这些判断环节让这个方法在面对真实复杂结构时稳得多。我自己现在处理环境激励数据时默认管线就是信号预处理 → 过一遍SSI-COV扫阶次 → 画稳定图 → MAC交叉验证 → 再回头调整p和滤波参数复核一遍。这个流程走下来误报模态的概率很小。如果你刚开始接触这个方法建议就用文章里的两自由度仿真算例做起点先把 p、n_modes、数据长度对识别结果的影响逐一试一遍跑出感觉之后再上实测数据。仿真算例最大的价值是可以随时用理论值校验你的代码和参数配置——识别结果和理论对得上说明流程通了对不上优先检查Toeplitz的下标和SVD截断逻辑这两个地方是SSI-COV代码最容易埋bug的位置。这篇文章从方法选型讲到原理、代码、验证和踩坑覆盖了SSI-COV做多自由度系统模态参数识别的主要环节。如果你正在做一个环境激励下的结构响应分析项目这套思路应该可以直接迁移过去。遇到具体问题——比如数据太短、通道太少、模态分不干净——可以沿着稳定图这条路慢慢排查经验都是跑数据跑出来的多折腾几次比看多少篇文章都管用。