SSI-COV随机子空间模态识别：Matlab实现与工程实践

做结构振动测试的人，十有八九遇到过这种尴尬：测了一堆加速度响应数据，激励力却测不到。风、车流、环境微振动一直在激励结构，我们根本没有办法记录真实的输入。这时候想做模态参数识别，传统的频响函数法直接失效，只能用“仅输出”的方法。而在这一众方法里，基于SSI-COV（协方差驱动随机子空间识别）的多自由度系统模态参数识别，是我在Matlab里用得最顺手、也最愿意跟人推荐的一套。它不依赖人工激励，不用迭代求解，从响应数据的协方差统计量出发，走一遍SVD和特征值分解就能把模态频率、振型和阻尼比全部提出来。

这篇文章想把SSI-COV从原理到代码到工程注意事项完整讲透。内容主要面向做结构健康监测、振动测试、动力学仿真的研究生和工程师，同时也适合刚接触模态识别、想在Matlab里快速落地一套可用流程的初学者。文中会包含可直接运行的Matlab核心代码、参数选取经验、稳定图实现思路，以及一个三自由度系统的完整仿真验证算例。看完之后你可以直接拿这套东西去处理自己的实测数据，也可以根据自己的问题场景做二次开发。

1. 为什么单选SSI-COV：环境激励下的识别困局

1.1 没激励信号时，传统方法为什么玩不转

经典模态分析里最常用的方法，是同时采集输入力信号和结构响应信号，然后算频响函数（FRF），再从幅频曲线的峰值去对频率和阻尼。这个方法在实验室里非常成熟，激振器一装、力锤一敲，FRF干净漂亮，参数识别精度很高。但到了现场，问题就来了：风荷载多大？车流激励怎么测？环境微振动压根没有一个明确的激励源。你既不可能在桥墩上装力传感器，也不可能把整个结构停下来单独做一次受迫振动试验。

这种情况下，基于输入输出的方法全部失效。于是业内转向“仅输出”模态识别：我不需要知道激励是什么，只假设它是平稳随机信号，从结构响应的统计规律里把模态信息反推出来。

常用的仅输出方法有几类。频域分解法（FDD）思路简单，把响应功率谱密度矩阵做SVD，奇异值曲线的峰值对应模态频率，但它给出的阻尼比精度一般，只能作为初步参考。随机减量法（RDT）或NExT技术需要先从响应里提取自由衰减响应，再做时域拟合，可是提取过程中很依赖经验参数，稍不注意自由衰减信号就变形了。Eigensystem Realization Algorithm（ERA）本身是个好方法，但它需要脉冲响应或者自由衰减作为输入，等于说前面还得再接一道预处理。

相比之下，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或者两类方法对拍验证。

下面把几种常见仅输出方法的特性拉个表，方便对照选型：

2. SSI-COV的数理骨架：从协方差矩阵到系统矩阵

2.1 离散随机状态空间模型是整套识别的底座

SSI-COV的全部推导都建立在线性时不变系统的离散随机状态空间模型上：

x_{k+1} = 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_{k+i} * y_k^T]

这个序列有一个非常好的性质：在白噪声激励假设下，它和系统矩阵之间存在确定性的关系，可以写成

R_i = C * A^{i-1} * G

其中G = E[x_{k+1} * y_k^T]是状态与输出的互协方差矩阵。这意味着协方差序列本身就完整携带了系统矩阵A的信息。

实际计算时，我们用时间平均代替统计平均：

R_i = (1 / (N - i)) * Σ_{k=1}^{N-i} y_{k+i} * y_k^T

接着把这些协方差矩阵块拼成一个分块Toeplitz矩阵：

T = [R_i, R_{i-1}, ..., R_1; R_{i+1}, 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^T

SVD的奇异值按从大到小排列，如果系统阶次确实是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(l+1: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就得收一收。一个常见的起点是i=100到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秒的数据。如果数据太短，协方差序列在长滞后段几乎全是噪声，低频模态和阻尼比识别会非常差。

预处理方面，三个动作绝对不能省：去均值、去趋势、带通滤波。去均值是因为绝大多数传感器都有直流偏移，不去掉的话协方差里会混入一个巨大的常数项；去趋势是去掉测试过程中由于温度漂移、仪器零漂产生的缓慢变化趋势；带通滤波则是把关心的频带之外的噪声滤掉，提高整体信噪比。

三个关键参数的工程经验值放一张表：

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)*l+1:k*l, :) = Y(:, k:k+j-1); Yf((k-1)*l+1:k*l, :) = Y(:, k+i:k+i+j-1); end % Toeplitz矩阵: 未来块与过去块的协方差 T = (Yf * Yp') / j;

这一段注意循环里的索引。Yp第k行块对应Y从第k列开始连续取j列，Yf第k行块对应Y从第k+i列开始取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(l+1: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识别结果与理论值对比

调用核心函数，设置行块数i=200，系统阶次范围从2到80做稳定图。取稳定点聚类后的中位数，得到的结果与理论值对比如下：

频率识别非常稳定，误差基本在0.5%以内，这个精度主要得益于数据长度足够长、信噪比合理。阻尼比误差比频率误差大不少，尤其第三阶模态，识别值比理论值高了50%。这就是SSI-COV的典型特征：阻尼比识别方差大、整体偏高，需要使用稳定图多次平均来压低随机误差。

6.3 噪声水平和数据长度对识别结果的影响

我做了一组对比实验，分别改变信噪比和数据长度，观察第一阶模态频率和阻尼比的识别稳定性：

结论很清楚：频率识别对噪声和数据长度都不太敏感，但阻尼比非常敏感。如果你的目标是拿阻尼比做结构损伤识别之类的高精度分析，要尽量保证高信噪比、长数据记录，并且多次识别取平均。

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的选择、数据预处理、稳定图参数这些环节都需要实际操作者自己把握。建议你拿到新数据时，先跑一遍完整流程，再把参数微调一遍，对比两次结果是否一致，这样能避开大多数隐蔽的坑。后续如果遇到更复杂的场景，比如非线性结构或者非平稳激励，可以在这个基础上扩展时变子空间方法，也可以和机器学习方法结合做自动模态拣选，这些都是值得继续往下走的方向。