企业官网建设流程全解析

1. 为什么选SSI-COV：模态识别里的“只听输出”这条路

1.1 模态参数到底是什么

搞结构动力学、机械设备故障诊断或建筑桥梁监测的朋友，对“模态参数”这四个字一定不陌生。多自由度系统的自由振动，本质上可以分解成若干阶模态的叠加：每一阶模态都有一个固有频率、一个阻尼比、一组振型位移。频率决定了结构在什么频段上容易发生共振，阻尼比决定了振动衰减的快慢，振型则告诉你结构在各个测点上“怎么动”。这三个量合在一起，差不多就是线性结构动态特性的全部家底。

工程上怎么用？举个例子，你在桥梁上布了一排传感器，测到环境振动后，第一件事就是识别前几阶模态频率。如果某一阶频率在一年内持续下降，往往提示结构刚度发生了变化，后面再配合振型变化，就能进一步判断损伤大概在哪个区段。再比如楼顶的舒适度评估，横向低阶模态的阻尼比是重要输入参数。可以说，模态参数识别是结构健康监测这条链路上最基础也最关键的一环。

但难点在于，实际结构的激励往往不是你想测就能测到的。桥面上跑的车、吹的风，这些激励的真实时程你拿不到，甚至激励方向都不唯一。于是“纯输出识别”成了唯一出口，SSI-COV就是这条路上用得最多、也最不容易翻车的方法之一。

1.2 环境激励下的识别困局

实验室做锤击法模态测试时，锤子的力信号是已知的，测点响应是输，动辄就能算出频响函数H(f)，再通过峰值拾取或曲线拟合把模态参数拿下来。这套“输入-输出”识别方法很成熟，但搬到现场就会碰壁：你不能给一座大桥设计一个可控的激振器，更不能让运营中的风机按照你设定的频谱来振。环境激励不可控、不可测，但你想识别结构动态特性，怎么办？

一种粗暴做法是所谓“峰值拾取法”：把响应信号做FFT，频谱上看到峰，就认为那是一个模态频率。这个方法简单，但缺陷也明显：频率分辨率受制于窗长，阻尼比基本估不准，两个靠得很近的模态还会糊成一团。紧接着出现的频域分解法(FDD)通过奇异值分解改善了振型提取，但本质上还是在频域里做窗口选择，依然存在分辨率和窗泄漏问题。

真正在时域里把问题讲清楚的，是随机子空间识别。SSI-COV正是其中协方差驱动的分支。它把环境激励建模成随机过程，不需要知道力的具体时程，直接从测点响应的协方差出发，利用线性系统理论中的“Toeplitz矩阵可分解为观测矩阵和可控性矩阵之积”这个代数性质，把系统的状态空间矩阵A估计出来。有了A，频率、阻尼、振型一套全出来了。整个过程是严密的代数推导，不是靠峰值猜，所以对密集模态和噪声的容忍度都高得多。

1.3 SSI-COV和其它方法怎么选

纯输出方法里，除了SSI-COV，还有SSI-DATA、ERA、FDD等。我自己实际用下来，会把它们按应用场景分开看待。

先看SSI-DATA。它直接用原始响应数据构Hankel矩阵，理论上信息更完整，在某些极端情况下精度比SSI-COV好，但对存储和计算量的要求高不少，而且对原始数据里那种强噪声偶尔会不太稳。SSI-COV先把响应转成协方差矩阵再做子空间分解，相当于先做了一次统计平均，把不相关噪声在很大程度上洗掉了，矩阵规模也因此比SSI-DATA小很多。对一个10通道、50万点的数据，SSI-COV在普通笔记本上几秒就能出一张完整稳定图，SSI-DATA通常要多花几倍时间。

FDD在频域里做SVD，操作快、参数少，很适合做快速筛查，但阻尼比精度差，频率分辨率也受窗长限制。ERA理论上需要脉冲响应，现场没有，除非你先做频域转换再反变换，绕一圈不如直接用SSI-COV。峰值拾取就更不用说了，连频率也经常看不准。

所以我的默认选择是：快速摸底用FDD，正式出模态参数用SSI-COV，关键数据再拿SSI-DATA交叉验证。这个组合在多个项目里都稳得很。

2. 原理拆解：协方差、Toeplitz矩阵和子空间

2.1 把系统写成状态空间

SSI-COV的起点，是假设结构可以用一个离散线性时不变随机状态空间模型描述：

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

其中x是n维内部状态向量，A是系统矩阵，C是观测矩阵，y是你测到的m通道响应，w和v分别代表过程噪声和测量噪声，都假设为零均值白噪声。环境激励下，那个测不到的广义力被统一吸收进w里，你只面对象征着随机激励的噪声项。

为什么这套变量设定能导出模态参数？因为A的特征值和特征向量，直接对应的正是结构在离散时间下的极点。每对物理模态都会在A中贡献一对共轭复数特征值，复频率的模对应圆频率，实部对应衰减率。C的作用则在于把状态空间里的内部自由度投影到传感器测点位置，C乘以特征向量就得到测点上的振型分量。

也就是说，只要能估计出A和C，模态频率、阻尼比、振型全部都能读出来。SSI-COV的核心工作，就是想办法从输出响应序列里把A和C“逼”出来。

2.2 从输出信号到协方差矩阵

既然只有测点响应y_k，第一步是构造延时协方差矩阵：

R_τ = E[y_{k+τ} y_k^T]

实际操作中，用样本平均估计：

R_τ ≈ (1/(N-τ)) Σ_{k=1}^{N-τ} y_{k+τ} y_k^T

为什么协方差这么好用？关键在于白噪声的性质：白噪声在不同时刻是完全不相关的，所以它在τ≠0的延时协方差里几乎不贡献任何成分，只有τ=0时才有方差贡献。换句话说，你在计算延时协方差时，等于做了一次天然的降噪处理——传感器测量噪声和环境输入噪声里与系统响应无关的部分，在非零延时域被自动排除了。这一点是SSI-COV在抗噪性上优于直接时域方法的原因之一。

把所有延时都算出来后，把它们按行排成块Toeplitz矩阵：

T_{1|i} = [R_i, R_{i-1}, ..., R_1; R_{i+1}, R_i, ..., R_2; ... R_{2i-1}, R_{2i-2}, ..., R_i]

这里的i叫块行数，决定了用多少阶延时信息。这个矩阵看起来规规矩矩，背后却藏着一个关键代数结构：它恰好等于观测矩阵O_i和另一个可控性型矩阵的乘积。这是随机子空间方法最核心的漂亮结论。

2.3 SVD、观测矩阵和状态矩阵A

对T_{1|i}做奇异值分解：

T = U S V^T

因为T能够分解成O_i乘以某个矩阵，所以T的列空间和O_i的列空间在本质上是一致的。你截取S中前n个大奇异值，对应系统真实阶次，舍掉后面的小奇异值，就能得到截断后的观测矩阵：

O = U(:,1:n) * sqrt(S(1:n,1:n))

观测矩阵O具有分块结构，第一块恰好就是C，而它的“上下移位关系”内藏着A。因为O的每一块都是上一个块乘A，所以构造O_up为去掉最后一块行、O_down为去掉第一块行，就有：

A = pinv(O_up) * O_down

这一步是整个计算的灵魂：看似只用了矩阵移位，实际上把状态转移矩阵A从观测数据里代数地解出来了。

2.4 从特征值到频率、阻尼、振型

拿到A后做特征值分解，设某个离散特征值为λ_d。因为A是离散状态矩阵，与其对应的连续时间极点λ_c满足：

λ_c = ln(λ_d) / Δt

有了连续极点，三个参数直接读出：

f = |λ_c| / (2π) ξ = -Re(λ_c) / |λ_c|

振型则从观测矩阵C和特征向量V_A一起来：

φ = C * V_A

这里有个细节要注意：λ_c通常以共轭复数对出现，一对就对应同一阶物理模态的正负频率项，识别时只用正频率项即可，否则统计结构时会重复计数。另外，输出只有条件识别出的振型是不带质量归一化的，也就是说只能看相对振型，不能像试验模态那样直接做振型缩放。大多数监测和损伤识别场景要的就是相对振型，所以这不是大问题。

3. Matlab完整实现：从仿真数据到稳定图

3.1 先造一组三自由度仿真数据

写算法之前，我强烈建议先造一组带真实答案的仿真数据，否则你根本不知道算法是对是错。下面这个三自由度链式系统是我常用的测试算例，阻尼用比例阻尼模拟。

%% 三自由度弹簧-质量系统仿真 dt = 0.005; % 采样间隔 200 Hz t = (0:dt:60)'; % 60 秒 M = diag([2 1 1]); K = 1000*[1 -1 0; -1 2 -1; 0 -1 1]; alpha = 0.05; beta = 0.0002; C = alpha*M + beta*K; A_ss = [zeros(3,3) eye(3); -M\K -M\C]; B_ss = [zeros(3,3); M\eye(3)]; C_ss = [eye(3) zeros(3,3)]; D_ss = zeros(3,3); sys = ss(A_ss, B_ss, C_ss, D_ss); u = randn(3, numel(t)); % 三通道不相关的白噪声激励 y = lsim(sys, u, t)';

这里用连续时间状态空间生成响应，再用lsim做仿真，既保证了模态成分真实存在，又通过白噪声输入模拟了环境激励场景。你可以先用eig(K, M)算出理论频率，记下来，后面用来对比。

3.2 SSI-COV核心函数一步步写

下面是完整的核心函数，每一步都加了注释。需要说明的是，这个版本主要为了讲清楚原理，代码做了适度简化，工程上可以再优化。

function [fn, xi, phi, A, Cm] = ssicov(y, dt, i, n) % SSI-COV 协方差驱动随机子空间识别 % y : Ntime x Nchannels 输出信号 % dt: 采样间隔 (s) % i : 块行数 % n : 系统阶次 [N, m] = size(y); y = y - mean(y,1); % 1. 计算延时协方差 R_0 ~ R_{2i-2} R = zeros(m, m, 2*i-1); for k = 0:2*i-2 y1 = y(k+1:N, :); y2 = y(1:N-k, :); R(:,:,k+1) = (y1.' * y2) / (N-k); end % 2. 构造块 Toeplitz 矩阵 T = zeros(i*m, i*m); for p = 1:i for q = 1:i lag = i + p - q; % 对应 R_{i+p-q} if lag >= 1 && lag <= 2*i-1 T((p-1)*m+1:p*m, (q-1)*m+1:q*m) = R(:,:,lag); end end end % 3. 截断 SVD [U, S, ~] = svd(T, 'econ'); if n > size(S,1), n = size(S,1); end n = max(n, 2); % 至少保留 2 维，否则后面没法分解 Un = U(:,1:n); Sn = S(1:n,1:n); O = Un * sqrt(Sn); % 4. 利用移位性质计算 A 和 C O_up = O(1:(i-1)*m, :); O_down = O(m+1:i*m, :); A = pinv(O_up) * O_down; Cm = O(1:m, :); % 5. 特征值分解 + 离散转连续 [V_A, D_A] = eig(A); lam_d = diag(D_A); lam_c = log(lam_d) / dt; % 6. 提取模态参数 fn = abs(lam_c) / (2*pi); xi = -real(lam_c) ./ abs(lam_c); phi = Cm * V_A; % 7. 只保留正频率，并按频率升序排列 idx = imag(lam_c) > 1e-6; fn = fn(idx); xi = xi(idx); phi = phi(:, idx); [fn, order] = sort(fn); xi = xi(order); phi = phi(:, order); end

几点实现说明：

• 协方差计算用矩阵乘法y1.' * y2，比for循环逐点累加快得多，数据量大时尤其明显。
• Toeplitz 矩阵的索引是i + p - q，第一次写的时候我也经常搞反，建议你对照公式手推一遍第一行的R_i, R_{i-1}, ..., R_1。
• eig(A)默认会把小特征值清零，为了稳定我习惯用eig(A, 'nobalance')，对含噪声的A矩阵更靠谱。
• 阶次n不能太小，至少要比你预计的物理模态数乘2大一点，否则真实模态会被截掉。

3.3 稳定图怎么画

真实数据里，物理模态和数值伪模态混在一起，光靠单次n的选择很难判断谁真谁假。工程上的标准做法是画“稳定图”：让系统阶次从2逐步增加到40甚至更大，每个阶次都识别一遍，把极点画在“频率-阶次”坐标平面上。物理模态因为对模型阶次不敏感，会在对应频率位置形成一条整齐的垂直稳定带，而伪模态则东一个西一个，前后对不上。

figure; hold on; for n = 2:2:40 [fn_cur, ~, ~] = ssicov(y, dt, i, n); scatter(fn_cur, n*ones(size(fn_cur)), 10, 'filled'); end xlabel('频率 (Hz)'); ylabel('系统阶次 n');

只看散点还不够，更严格的稳定判据通常是三条：频率变化小于1%、阻尼比变化小于5%（或10%）、MAC值大于0.95。我习惯在代码里加一个循环，把连续三阶以上都满足判据的极点标记为“稳定点”，画成实心大圆点。稳定点连成的竖线对应的频率，十有八九就是真实模态。

3.4 识别结果和理论值对比

仿真数据的好处是有标准答案。用eig(K, M)算出理论圆频率后，除以2π得到理论固有频率。我的经验是：在无噪声或低噪声比例阻尼仿真里，SSI-COV识别出的频率误差能控制在0.1%以内，阻尼比误差在5%~10%左右，振型用MAC值衡量，通常都在0.99以上。

如果你发现频率对，但振型MAC值不高，先检查通道顺序和C矩阵的行对应关系——传感器编号顺序错了，振型必然对不上。这个坑我踩过不止一次。

4. 参数选择与避坑：这些坑我替你踩过

4.1 块行数和模型阶次真的不能乱拍

SSI-COV有两个要手动定的关键参数：块行数i和模型阶次n。

块行数i决定Toeplitz矩阵里用了多少延时协方差信息。i太小，系统动态信息装不下，尤其低频模态的长周期相关抓不住；i太大，则会把噪声的统计残留也装进去，奇异值间的界限变得模糊，伪模态增多，计算量也变大。一个常见经验是让最大延时(2i-1)*dt覆盖最低阶感兴趣模态周期的2到4倍。还有一个更直接的现场做法：画一下第一通道的自相关函数，看它衰减到噪声底的位置大致对应多少秒，这个时间长度对应的延时数乘m，就可以作为i的参考量级。我通常会多试几个i值，稳定图上模态带最清晰的那个i就是好选择。

模型阶次n的选取则要靠“过参数化”思维。真实物理模态每阶占两个极点（共轭对），假如系统有3阶真实模态，n至少取6。但现场数据里总有噪声和未建模的动态，所以要把n设得比实际需要的更大，比如n从2一直扫到40。真实模态在稳定图上会形成稳定竖线，伪模态则散乱分布。千万别指望一次就能把n选准，稳定图存在的意义就是让你不需要事前猜准n。

4.2 数据预处理比算法本身更影响结果

很多新手把精力全放在算法和代码上，却忽视了数据本身。SSI-COV假设激励是平稳随机过程，站点数据里普遍存在的趋势项、传感器漂移、工频干扰都会破坏这个假设。

第一步永远是去除均值和线性趋势，否则协方差矩阵会被趋势项主导。第二步是检查频谱，看关心模态频段内信噪比是否充足。如果高频噪声太大，先做低通滤波再重采样；如果只关心某一段频率的模态，可以用带通滤波把频段切出来，但要记得用filtfilt做零相位滤波，普通IIR滤波器会带来相位畸变，把振型搞坏。滤波器阶数不必太高，过度滤波反而会把真实模态削平。

第三个容易被忽略的问题是数据的平稳性。现场如果有突然的重车经过，响应瞬时非平稳，直接整段识别容易出错。我的做法是先把数据切成1到2分钟的分段，逐段独立识别，再看各段结果的分散度。稳定段之间频率结果如果非常接近，就取中位数作为最终估计；如果某段结果和其他段差得离谱，直接扔掉。

4.3 阻尼比为什么总是最不可靠的一个

频率和振型在SSI-COV里通常能识别得很漂亮，唯独阻尼比总让人捏把汗。原因在数学上很清楚：阻尼比由极点实部决定，而低频模态的实部相对虚部来说往往只有百分之几甚至千分之几，任何一点噪声扰动都会让这个微小量产生大比例误差。我在实测数据里见过同一段数据切两半，前半识别阻尼比2.1%，后半识别3.6%，频率倒是稳如老狗地停留在同一根竖线上。

所以我对阻尼比的建议是：不要指望单次识别定终身。做多段数据独立识别，画阻尼比直方图，看它中心落在哪、离散度有多大。如果直方图明显双峰，通常说明数据里混了几种不同的激励状态，要回到时域去查原因。另外，凡是识别出负阻尼比或阻尼比大于10%的极点，先别急着当成异常结构特性——大概率是伪模态或阶次选择不当，别把噪声的锅甩给结构。

5. 常见问题排查实录

5.1 快速排查表

这张表不是从哪本书上抄的，是我把那些年调试SSI-COV遇到的典型问题攒出来的，每行背后都对应过一次真实的折腾。

5.2 两个现场案例复盘

第一个案例是某实验室里的一根简支钢梁试验。6个测点、采样500Hz、采集了20分钟环境激励。第一版稳定图出来，22.5Hz附近有一条漂亮稳定线，但旁边25Hz处也有一条若隐若现的短线。一开始我以为是结构第二阶模态，后来把频段放大、增大块行数，发现那条短线其实来自实验室风扇的周期干扰，并不是结构模态。处理方法是把识别频段上限改到30Hz以下，并检查相干函数确认信号的物理来源。所以稳定图上任何可疑竖线，都要回到频谱和时域去交叉验证，不能只信算法。

第二个案例是某多层建筑的环境振动测量。九个通道同步采集，前几次识别频率很稳定，但阻尼比散得离谱，从0.8%到3.2%都有。后来我把数据按交通高峰和非高峰时段分开，发现非高峰段的阻尼比离散度明显小，而且中位数更接近预期值。这说明数据段里混入了大幅值的非平稳扰动。最后我采取分段识别+中位数稳健估计的方案，整套结果的可用性立刻上来了。

6. 最后聊几个实用细节

如果你准备把SSI-COV用到真实数据上，我最后再补三条个人经验。

第一条，算法跑通后先别急着换复杂工地。改参数时一次只动一个变量，比如先固定n扫i，再固定i扫n，看稳定图的变化规律。这样你能形成对算法的直觉，分得清哪些变化是参数引起的，哪些是数据里的真实特征。

第二条，保存SVD的奇异值序列。截断之前先画一下奇异值分布，如果前几个奇异值和后面之间有明显“台阶”，说明系统的有效阶次基本在那里，这个信息能帮你划定n的扫描区间，省掉很多盲目尝试。当然现场数据噪声大，“台阶”不一定清晰，但至少是个有价值的线索。

第三条，也是我踩过很多次坑后的体会：SSI-COV不是一个“数据一扔、参数一出”的黑箱工具。你花在数据预处理、稳定图判读、多段交叉验证上的时间和精力，往往比算法本身的实现要多得多。但正是因为这些前置工作做到位了，它才能在你的现场数据里给出真正可靠的结果。

做模态识别就像看病，逆问题是第一步，确认这步可靠，后面的诊断、预警、决策才有意义。希望这篇文章能帮你把SSI-COV这个工具箱真正用顺，让结构“开口说话”。