企业官网建设流程全解析

做结构动力学的人，迟早会遇到一个绕不开的问题：结构上明明测到了振动，却不知道是谁在敲它。桥梁上的车流、风机塔筒的风载、楼板上的人群，这些激励在测试时都测不到，可模态参数还是得认认真真识别出来。SSI-COV（协方差驱动的随机子空间识别）就是在这种“输出only”条件下干活的老牌方法：没有力信号，只用加速度响应，一次性把多自由度系统的模态频率、振型和阻尼比全提出来。我这两年把整套流程用Matlab跑过很多遍，从仿真模型到实测数据都踩过不少坑，这篇就按我的实操路径来写：为什么选它、原理怎么落地、代码怎么写、坑在哪里。适合刚接触运行模态分析、或者正准备用Matlab实现SSI-COV的朋友参考。

1. 为什么选SSI-COV：环境激励下的方法选型

1.1 识别困局：激励测不到，只能靠响应

传统试验模态分析的路子是锤击法或者激振器扫频，激励信号可控可测，频响函数一算，模态参数自然就出来了。但现场结构不是实验室，桥梁的风和车流、风机的随机风载、厂房的人行载荷，这些激励既无法精确测量，也没法重复，连参考通道都没有。

这时候只能走“输出only”路线。常见思路分两类：一类在频域做，从响应功率谱里挑峰值，代表作是峰值拾取法和频域分解（FDD）；另一类在时域做，用响应时间序列反推系统状态方程，代表作是ERA、ARMA以及随机子空间识别家族。时域方法不用做谱估计，规避了窗函数和频谱泄漏这一堆麻烦，而且能直接给出阻尼比，这正是频域方法最薄弱的地方。SSI-COV又是时域方法里概念最直观、计算量最可控的，所以不少Modal软件的内核都愿意用它。

1.2 SSI-COV和SSI-DATA、FDD的取舍

随机子空间识别其实有两条路线：数据驱动的SSI-DATA和协方差驱动的SSI-COV。两者共用同一套状态空间理论，差别在于怎么从响应里把系统矩阵抠出来。SSI-DATA直接对Hankel矩阵做QR分解再投影，数学上更稳，阻尼估计偏差也小一些；SSI-COV先把数据压成协方差Toeplitz矩阵再做SVD，计算量小、内存省，但协方差尾部误差会直接影响阻尼精度。

实际工程里怎么选？如果数据平稳、长度够，SSI-COV完全够用；如果信号噪声大、数据又紧张，我优先考虑SSI-DATA。FDD更适合快速摸谱峰，阻尼比只能靠半功率带宽凑合，密集模态下还容易糊。我最常用的组合是：FDD快速过一遍谱峰，SSI-COV精算频率和振型，阻尼比作为参考值。因为系统矩阵A的阶数直接对应模态数量，一次特征值分解把所有模态拿全，这就是SSI-COV对多自由度系统最友好的地方。

注意：SSI类方法的前提假设是激励近似为白噪声、系统线性时不变、输出平稳。现场激励如果偏离这个假设太多，任何输出only方法的结果都要打折。

2. SSI-COV的原理拆解：从响应协方差到系统矩阵

这一节把数学拆成五步走。别被SVD和状态空间吓到，本质上都是线性代数的常规操作。

2.1 多自由度系统的状态空间表达

多自由度系统的运动方程是 Mq'' + Cq' + Kq = f(t)。二阶方程之间是耦合的，为了用线性系统理论，把位移和速度拼成状态向量 x = [q; q']，写成状态空间形式：

x' = A_c x + B_c f，y = C x

采样后离散化成：

x_{k+1} = A x_k + w_k，y_k = C x_k + v_k

这里A是我们最关心的系统矩阵，频率、阻尼、振型全都藏在它的特征值和特征向量里。SSI-COV的目标，就是从输出数据y_k里估计出A和C，完全不需要知道输入f。你可以把A理解成一个“动力学编码器”，记录了系统每步如何演化，我们只能看到它的输出屏幕，但可以通过输出的统计关联把这个编码器反推出来。

2.2 第一步：用Hankel矩阵组织数据

识别前要把原始响应重新排列成特殊形式的矩阵——Hankel矩阵。以两倍块行数切开，上半块是“过去”输出，下半块是“将来”输出，每个块里都放着所有测点的时间序列。

为什么要切成过去和将来？因为系统的动态特征会体现在“过去的输出和将来的输出之间的相关性”里。块行数取i，那么Hankel矩阵维度是2il行×（N-2i+1）列，l是测点数，N是每个测点的样本数。过去块Y_past和将来块Y_future都是i*l行。这一步开始前去掉均值和线性趋势非常关键，否则直流分量的能量会污染整个协方差矩阵。

2.3 第二步：由协方差生成Toeplitz矩阵

随机状态下，输出的协方差序列满足 R_m = E[y_{k+m} y_k^T] = C A^{m-1} G，其中G是状态和输出之间的特定统计矩。换句话说，协方差序列和系统矩阵A之间存在确定的映射关系。

将过去块和将来块做一次外积平均，就得到Toeplitz矩阵 T = Y_future * Y_past' / j。T是分块Toeplitz矩阵，每个块本质上是不同滞后的协方差。由于A的所有极点信息都会“藏”在这种延迟相关结构里，T相当于把系统的特征在统计层面集中了起来。计算T有两种写法：一种是循环算R_0, R_1…再填块，一种是用Hankel外积法。后者代码简单，不容易把块填错，实测两者结果一致。

2.4 第三步：SVD分解提取可观测矩阵

T满足结构T = Γ_i · O_i，其中Γ_i是可观测矩阵，由C、CA、CA²…堆叠而成；O_i是随机可控制矩阵。对T做奇异值分解 T = U S V'，取前n个最大奇异值截断，令Γ_i = U_1·sqrt(S_1)，就把T的左右因子之一解出来了。

奇异值的大小体现了对应模态在响应中的能量占比，所以奇异值谱本身有一定参考意义。但现实中受噪声和模态能量差异影响，靠奇异值“断裂处”定阶并不可靠，更稳妥的还是稳定图。这里还要注意：n是模型阶数，它等于实际模态阶数的2倍，因为每阶实数模态对应一对共轭极点，所以参数扫描时要按2递增。

2.5 第四步：最小二乘反推系统矩阵A

可观测矩阵有一个重要性质：Γ_i去掉最上面一块得到Γ_up，去掉最下面一块得到Γ_down，两者满足Γ_down = Γ_up · A。直接两边左乘pinv(Γ_up)，就能得到A的最小二乘估计。

这一步原理不复杂，但数值上要小心。当块行数i选得太小，Γ_up本身病态时，pinv出来的A可能带明显误差。实测经验是：把i取得足够大，让Γ_up有足够的行来约束解，比强行增加模型阶数n更有用。

2.6 第五步：特征值分解得到频率、阻尼比和振型

最后把A做特征值分解：A ψ = λ ψ。离散特征值λ要映射回连续时间系统，用 s = log(λ)/Δt。连续极点s的实部和虚部分别对应衰减率和振荡频率：

• 模态频率 f = |s| / (2π)
• 阻尼比 ζ = -Re(s) / |s|

注意这里算出来的频率本质上是无阻尼固有频率，低阻尼时和阻尼频率几乎没差别，不用纠结。振型提取要做两步：A的特征向量反映的是状态空间里的模态向量，要乘以输出矩阵C才转到物理测点自由度上，也就是φ = C·ψ。比例阻尼系统里φ本来是复向量，但各元素相位接近0或π，工程上取实部再归一化，就是常见的位移振型。

3. Matlab全流程实现：2自由度仿真模型跑通SSI-COV

3.1 仿真系统的搭建：质量、刚度与阻尼

先用一个2自由度弹簧-质量系统把流程走通。参数我选得比较常规：m1=8kg，m2=5kg，三个弹簧刚度分别为k1=2000、k2=2800、k3=3000 N/m。这样算出来的两阶固有频率大约在3Hz和6Hz，分离得比较开，适合先验证算法流程。

%% 1. 定义2自由度系统参数 m1 = 8; m2 = 5; k1 = 2000; k2 = 2800; k3 = 3000; M = [m1 0; 0 m2]; K = [k1+k2, -k2; -k2, k2+k3]; % 理论模态：求解广义特征值问题 [Vth, Dth] = eig(K, M); w_th = sqrt(diag(Dth)); f_th = w_th / (2*pi); % 瑞利阻尼：指定前两阶阻尼比 zeta1 = 0.02; zeta2 = 0.03; w1 = w_th(1); w2 = w_th(2); beta = 2*(w2*zeta2 - w1*zeta1) / (w2^2 - w1^2); alpha = 2*w1*zeta1 - beta*w1^2; Cd = alpha*M + beta*K;

瑞利阻尼里阻尼比和频率的关系是 ζ_i = α/(2ω_i) + βω_i/2，两个方程联立就能解出α和β。我见过不少初学者直接乱填Cd矩阵，结果仿真出来的衰减特性根本不对，后面识别出的阻尼比自然对不上。先用两阶指定阻尼比反推α、β，是最可控的做法。

3.2 白噪声激励与响应生成

激励用两个自由度上的高斯白噪声力信号，按100Hz采样，记录200秒。这个设置模拟环境激励的宽频特性，也顺便考验算法从噪声里挑出模态的能力。

%% 2. 生成白噪声激励下的响应 fs = 100; dt = 1/fs; T = 200; t = (0:dt:T).'; N = length(t); n = 2; A_ss = [zeros(n), eye(n); -M\K, -M\Cd]; B_ss = [zeros(n); inv(M)]; C_obs = [1 0 0 0; % 观测位移，两个自由度 0 1 0 0]; D_obs = zeros(2,2); sys = ss(A_ss, B_ss, C_obs, D_obs); rng(2025); F = randn(N, 2); Y = lsim(sys, F, t); % N×2，列是测点 % 加一点点测量噪声，更接近实际 noise_std = 0.02 * mean(std(Y, 0, 1)); Y = Y + noise_std * randn(size(Y));

这里我故意用的是位移观测而不是加速度观测，因为位移的振型最直观。实际测试多用加速度计，这时观测矩阵要换成 C_obs = [-M\K, -M\Cd]，振型就是加速度振型，思路完全一样。

3.3 SSI-COV主函数：从响应矩阵到模态参数

下面这个函数就是SSI-COV的主干。入参是输出响应矩阵Y（l×N）、块行数i、模型阶数n、采样间隔dt，出参是频率、阻尼比、复振型以及中间量。

function [f, zeta, phi, A, lam, s_c] = ssi_cov(Y, i, n, dt) % SSI-COV 协方差驱动随机子空间识别 % Y: l×N，l为测点数，N为样本数 % i: 块行数 % n: 模型阶数（真实模态数的2倍左右） % dt: 采样时间间隔 [l, N] = size(Y); Y = Y - mean(Y, 2); % 去均值，去除直流分量 if N < 2*i*l + 1 error('数据长度不够，请减小块行数i'); end % 构建Hankel矩阵，切成过去块和将来块 r = i*l; H = zeros(2*r, N - 2*i + 1); for k = 1:2*i H((k-1)*l+1:k*l, :) = Y(:, k : k + (N - 2*i)); end Y_past = H(1:r, :); Y_future = H(r+1:end, :); % Toeplitz矩阵 = 协方差矩阵的块状堆叠 T = (Y_future * Y_past') / (N - 2*i + 1); % 奇异值分解，取前n个奇异值截断 [U, S, ~] = svd(T, 'econ'); S = diag(S); U1 = U(:, 1:n); S1 = diag(S(1:n)); Gam = U1 * sqrt(S1); % 由可观测矩阵的分块性质反推系统矩阵A l2 = l*(i-1); Gam_up = Gam(1:l2, :); Gam_down = Gam(l+1:end, :); A = pinv(Gam_up) * Gam_down; % 特征值分解 [Psi, D] = eig(A); lam = diag(D); % 连续极点 s_c = log(lam) / dt; f = abs(s_c) / (2*pi); zeta = -real(s_c) ./ abs(s_c); % 振型：状态空间特征向量经过观测矩阵转到物理自由度 Cmat = Gam(1:l, :); phi = Cmat * Psi; end

几个容易忽略的细节：

• 为什么T用Y_future * Y_past'而不是反过来？因为按这个方向构造的Toeplitz块是输出协方差，对应Γ·O的结构，SVD之后可以直接拆出可观测矩阵。反过来的话结构关系就乱了。
• pinv(Gam_up) 用的是伪逆，因为Gam_up不是方阵，只有在最小二乘意义下能解出A。
• log(λ)这一步很多人漏了，直接把离散特征值当连续极点用，频率和阻尼会完全跑偏。
• 特征值分解后得到的Φ，行数等于状态维数n。但物理测点数l一般小于n，所以要先乘以Cmat，才能得到在每个传感器位置上的振型幅值。

3.4 筛选振动极点与结果对比

n如果取10，A就是10×10，特征值会出来10个极点，但真实系统只有4个有物理意义的极点（两阶模态乘以共轭对），其余都是噪声模式。筛选时把非振动极点剔除：

% 调用SSI-COV n_order = 10; [f, zeta, phi, A, lam, s_c] = ssi_cov(Y.', i, n_order, dt); % 筛选：保留虚部不为0、频率在合理范围内的极点 vib = abs(imag(s_c)) > 1e-3 & f > 0.1 & f < 0.9*(fs/2); f = f(vib); zeta = zeta(vib); phi = phi(:, vib); % 按频率排序 [~, idx] = sort(f); f = f(idx); zeta = zeta(idx); phi = phi(:, idx);

用固定随机种子跑完一组，我得到的典型结果和理论值对比如下：

MAC（Modal Assurance Criterion）是振型相关性的标准判据，接近1说明识别振型和理论振型高度一致。频率识别得很准，阻尼比有一定偏差，这完全是SSI-COV的正常水平。不要指望它把2%的阻尼识别成精确的2.000%，能到1.8%~2.2%这个区间已经可以接受。

3.5 稳定图定阶的思路

前面几次调用ssi_cov时，n是拍脑袋定的。实际工程中不可能预先知道模态阶数，所以要把n从2往上扫，把每个阶次识别出的频率画成散点图，能叠加成一条“竖直轨道”的，才是真实模态。这就是稳定图。

nmax = 20; figure; hold on; for n = 2:2:nmax [f_t, zeta_t, phi_t, A_t, lam_t, s_t] = ssi_cov(Y.', i, n, dt); keep = abs(imag(s_t)) > 1e-3 & f_t > 0.1 & f_t < 0.9*(fs/2); scatter(f_t(keep), n*ones(sum(keep),1), 8, 'filled'); end xlabel('Frequency [Hz]'); ylabel('Model order n');

画出来之后，你会在3Hz和6Hz附近看到两条笔直的稳定线，其他位置是散乱的噪声点。稳定线的判定一般看三个指标：频率相对变化小于1%、阻尼比绝对变化小于0.005、MAC大于0.95。能同时满足这三条的极点数最多、连续阶次最多的位置，就是推荐的定阶位置。稳定图看起来像个二维直方图，这也是整个流程里最“见真章”的一步。

4. SSI-COV能用在哪儿：应用场景与组合打法

4.1 最适合SSI-COV的场景条件

SSI-COV不是什么都能干，但它适合的场景在工程里相当常见：

• 结构长期处于平稳环境激励下，比如桥梁正常通行、风机稳定运行、楼板正常使用。
• 关心的模态集中在低频段，比如大跨结构的0.2~5Hz频段，采样需求不高。
• 传感器通道数不算太多，几十通道以内SSI-COV的内存压力可以忽略。
• 需要批量处理长期监测数据，SSI-COV的流程固定，方便自动化。

反过来，如果激励非常不平稳，比如强冲击占主导，或者信号里含有明显的窄带谐波分量，SSI-COV就很容易把那些谐波当成模态识别出来，稳定图上出现“假轨道”。这时候先得确认数据里哪些是物理模态，哪些是激励特征。

4.2 与FDD、SSI-DATA、贝叶斯方法的组合

我自己做项目时不会只依赖一种方法，SSI-COV更多是承担“快而全”的角色。各类方法的定位如下：

工程打法通常是：FDD先看谱，圈出明显的谱峰；SSI-COV扫稳定图，把候选模态锁死；遇到阻尼比有争议的关键模态，再用SSI-DATA或者贝叶斯方法出正式报告值。这样每一步都有交叉验证，不容易被单一算法的缺陷带偏。

4.3 从离线分析到在线自动监测

SSI-COV的固定流程很适合嵌入到结构健康监测系统里。做法是：数据流按固定时长滑窗，每个窗口跑一次SSI-COV，自动画稳定图、聚类稳定极点，再把识别的频率和阻尼比推送到趋势曲线。阻尼比趋势对结构损伤比较敏感，所以很多桥梁和风力机的监测系统都用这种“自动随机子空间”套路。

自动化要注意两点。第一，滑窗长度不能太短，一般至少要有几十个结构自振周期，否则阻尼比估计方差很大。第二，稳定图的自动判定阈值要保守一点，宁可少报真模态，也别多报假模态，因为监测系统误报警比漏报警更麻烦。

5. 实操避坑实录：阻尼、稳定图与数据处理的常见坑

5.1 阻尼比为什么总不像直觉想的那么准

SSI-COV对频率和振型的估计通常很好，但阻尼比天生就容易波动。原因在于协方差的滞后越长，样本方差越大，而Toeplitz矩阵远处那些块恰恰对阻尼估计影响很大。低阻尼结构尤其明显：衰减慢，协方差尾部长，误差被放大。

对策方面，一是数据长度要够，别拿几秒钟的数据去估0.5%的阻尼；二是块行数i别取太大，否则远处噪声块白白拉高方差；三是多跑几次稳定图，看阻尼比在哪个区间内稳定，报一个范围而不是单点值。

5.2 稳定图上凭空多出“稳定轨道”

稳定图上有稳定线不等于一定是物理模态。信号里的正弦干扰、传感器局部共振、数据趋势没有被完全去掉，都可能在特定频率上形成假稳定线。判断方法有三招：

• 把响应功率谱叠在稳定图旁边，物理模态基本都对应谱峰，而假轨道不一定有对应的谱峰。
• 看看该极点的复平面分布，阻尼比为正且落在大致合理的区间才可信。
• 改变数据段重算，物理模态应该稳定出现，假轨道会在不同数据段里漂移。

5.3 块行数i和模型阶数n怎么选

这是SSI-COV被问得最多的问题。我个人的经验公式：i取感兴趣的最高阶模态数的2~5倍，再乘个安全系数。比如关心的模态不超过5阶，i取20到40之间比较合适。i太小，可观测矩阵的行数不够，A矩阵解不稳；i太大，Toeplitz矩阵里的远滞后协方差噪声增多，阻尼估计恶化。

n的扫描范围从2开始，步长2递增，最大值取期望模态数的4~6倍。这样做既不会漏掉真实模态，又能给噪声模式留出空间，稳定图才看得出哪些是“站得住”的极点。

5.4 数据预处理的三板斧与第四板斧

去均值、去趋势、异常段剔除是基本操作。特别强调第四板斧：带通滤波。如果你的分析频带只有0.5~10Hz，却拿原始信号直接算，那么工频干扰、高频仪器噪声、低频漂移都会优先占据Toeplitz矩阵的主导奇异值，把真实模态挤到后面去。滤波不会改变真实模态的位置，但会显著提升信噪比。

滤波时要小心相位失真，最好用零相位滤波（filtfilt），这样不会因为相位移动影响振型的相对关系。

5.5 采样率、单位、数值病态这些细节

采样率不够，模态频率靠近或超过Nyquist频率，特征值的主辐角会接近±π，log映射之后频率就错位了。经验法则是采样率比最高目标频率高5到8倍以上，宁可多采也别少采。如果只关心低频，先抗混叠滤波再降采样，SSI-COV算起来更快，稳定图也更干净。

单位问题也很隐蔽。位移用米还是毫米，响应矩阵的数值尺度差1000倍，SVD的奇异值会差出好几个数量级。别让数据尺度过分悬殊，必要时统一单位，否则截断判别和pinv的数值条件全会受拖累。

5.6 避坑速查表

跑过几轮从仿真到实测的完整流程之后，我的体会是：SSI-COV不是“一个脚本打天下”的工具，它对数据的脾气很敏感。但只要你把平稳段挑好、滤波调对、稳定图反复核对，它给出的频率和振型通常非常可信。阻尼比就别太较真，能到工程判断的精度已经值得庆幸。最后分享一个小技巧：每次跑参数扫描之前把随机种子固定住，这样不同阶次的识别结果可以在同一份数据上互相对比，稳定图看起来才不散乱，代码也方便复现。