企业官网建设流程全解析

我前阵子帮一个做土木的朋友处理过一串实测楼板振动数据，慢慢意识到一个问题：很多搞结构的人手里晃着好数据，却对"怎么把模态参数干净地摘出来"这一步犯怵。频域峰值法简单但主观，拟合传递函数又对激励有硬性要求，到了现场实测这种环境激励占主导的场景，传统方法常常翻车。这次就把我在多自由度系统上用SSI-COV（协方差驱动随机子空间识别）做模态参数识别的完整流程捋一遍，包括模态频率、阻尼比、振型的Matlab代码实现，以及我踩过的一些坑。

这套方法适合下面这些人看：要做运行模态分析（OMA）的工程研究人员、做结构健康监测的从业者，以及刚入门模态分析但不想只停留在"调用工具箱"层面的学生。全文用一个三自由度剪切型结构做例子，白噪声激励模拟环境振动，从仿真数据生成到SSI-COV原理，再到Matlab代码每一步的实现，最后用稳定图筛选物理模态，全程可复现。

1. 模态参数识别到底在解决什么问题

先搞清楚一个最基础的问题：我们为什么要做模态参数识别。一个多自由度系统的运动方程可以写为：

Mẍ(t) + Cẋ(t) + Kx(t) = f(t)

M、C、K分别是质量、阻尼、刚度矩阵。直接去测M和K是不现实的，结构太大、边界条件太复杂，有限元模型也总有误差。但系统在振动中会自然暴露出它的固有属性——模态频率、阻尼比、振型，这三个参数统称模态参数。它们本质上是由M、C、K共同决定的特征信息，反过来通过实测响应把这三个参数提取出来，就是模态参数识别做的事。

传统方法里，锤击法或激振器法需要已知激励力，测力信号和响应信号做频响函数估计。这在实验室里没问题，到了实际工程现场，你很难对一座桥或一栋楼施加可控的激励，更多的场景是利用环境激励——风、地脉动、交通荷载。此时激励不可测，只能依靠响应数据做识别，这催生了运行模态分析（Operational Modal Analysis, OMA）这一大类方法。SSI-COV就是OMA领域里公认精度高、抗噪能力强的代表方法之一。

还有一点值得说清楚：模态参数识别不只是为了"测出来好玩"。后续的有限元模型修正、损伤识别、结构健康监测，甚至舒适度评价，全部建立在这三个参数的准确性之上。阻尼比尤其敏感，它数值小、对噪声敏感，很多近似方法算出来都是"能看但没法用"。SSI-COV在阻尼比识别上的表现，是它在工程上受欢迎的重要原因。

2. 为什么是SSI-COV：方法比较与原理拆解

2.1 运行模态分析的几类主流思路

运行模态分析在工程上的方法基本分三大派系：

• 频域方法：以峰值拾取法（Peak Picking）、频域分解法（FDD）为代表。FDD通过对响应功率谱密度矩阵做SVD分解，用奇异值曲线峰值定位频率。优点是快、直观，缺点是阻尼比估计粗糙，频率分辨率受FFT影响，密集模态容易混在一起。
• 时域方法：以ITD、STD、特征系统实现算法（ERA）、随机子空间识别（SSI）为代表。直接从时域响应中提取状态空间模型，再计算模态参数。精度高但计算量相对大。
• 时频方法：小波变换、HHT等，适合时变系统，这里不多说。

SSI又分为数据驱动的SSI-DATA和协方差驱动的SSI-COV。SSI-DATA直接对时间序列数据做QR分解和SVD，数值稳定性好；SSI-COV则先从响应数据构造协方差矩阵，再对这个矩阵做SVD分解。两者本质上是等价的，SSI-COV的计算量更小，代码实现也更直观，这正是我在Matlab里优先选择它的原因。

2.2 SSI-COV的数学原理：从响应数据还原状态空间模型

SSI-COV的核心思想绕不开一个关键前提：一个线性时不变系统，在环境激励下，其输出响应可以用离散状态空间模型描述：

x(k+1) = A·x(k) + w(k) y(k) = C·x(k) + v(k)

其中x(k)是系统内部状态向量，y(k)是测量输出向量，w(k)和v(k)分别是过程噪声和测量噪声，假设都是零均值白噪声。A是系统矩阵，它的特征值里藏着系统的频率和阻尼；C是输出矩阵，它和A的特征向量结合就能还原振型。

所以整件事变成了：在激励不可测的情况下，仅凭y(k)把A和C估计出来。SSI-COV的巧妙之处在于它绕开了不可测输入，转而利用输出数据的协方差信息来构造系统矩阵的估计。

推导路径是这样的：

第一步，构建Hankel矩阵并划分过去和未来块。把通道数为L、时长为N的响应数据y(k)排列成一个分块Hankel矩阵，上半部分叫"过去"Yp，下半部分叫"未来"Yf：

Yp = [y(0) y(1) ... y(N-2i)] [y(1) y(2) ... y(N-2i+1)] ... [y(i-1) y(i) ... y(N-i-1)]

Yf = [y(i) y(i+1) ... y(N-i)] [y(i+1) y(i+2)... y(N-i+1)] ... [y(2i-1) y(2i) ... y(N-1)]

这里的i是块数，每个块包含L行，过去和未来各i块。

第二步，用未来块和过去块的互协方差构造Toeplitz矩阵。数学上可以证明，输出协方差矩阵R(k) = E[y(k)·y(0)^T]与系统矩阵之间存在关系。把所有需要的协方差堆叠起来，得到Toeplitz矩阵T(1|i)，它的表达式恰好可以分解为可观测性矩阵O_i和可控性矩阵Γ_i的乘积：

T(1|i) = Yf · Yp^T / N = O_i · Γ_i

O_i矩阵里含C和A的幂次，这就是为什么后续能从它身上还原系统矩阵。

第三步，对Toeplitz矩阵做SVD分解。T = U·S·V^T，奇异值从大到小排列。理论上前2n个奇异值（n是模态阶数）对应真实模态，剩下的奇异值接近零，反映了噪声。取前2n个奇异值，把U、S截断为U1、S1，就得到了降秩后的Toeplitz矩阵。

第四步，从SVD结果还原C和A。可观测性矩阵O_i = U1·S1^(1/2)。输出矩阵C直接取O_i的前L行——因为可观测性矩阵的第一行块就是C。系统矩阵A利用可观测性矩阵的位移结构，用最小二乘求：

A = pinv(O_i(1:(i-1)L, :)) · O_i(L+1:iL, :)

第五步，对A做特征值分解，从特征值里解出频率和阻尼比。这一套流程逻辑非常严密，每一步都有扎实的线性代数支撑，这也是它比频域峰值法更让人放心的原因。

2.3 几个关键点的直观理解

这里用大白话解释一下为什么SVD能帮我们定阶。想象T矩阵的信息分为"真实振动"和"噪声"两部分，前者在奇异值谱上表现为较大的奇异值，后者散布在小奇异值上。做SVD相当于把矩阵按"信息重要程度"重新分解排序，只要真实模态的奇异值明显大于噪声奇异值，截断点就一目了然。实际数据如果信噪比不高，奇异值可能平滑衰减没有明显台阶，这时候就需要稳定图辅助判断（后面专门讲）。

另一个关键点是：每阶物理模态对应一对共轭复数极点，所以在状态空间里系统阶次是2n而不是n。代码里如果扫描到5阶模态，对应系统阶数应该设为10。新手经常在这里犯迷糊。

3. 多自由度系统仿真数据准备

3.1 三自由度剪切结构模型的建立

为了验证SSI-COV的实现效果，我建了一个经典的三层剪切型结构模型，模型简图就是三个质量块串联，层间刚度和阻尼集中。取质量m1=m2=m3=5000 kg，层间刚度k1=k2=k3=2×10^6 N/m。

质量矩阵和刚度矩阵在Matlab中这样构造：

m = 5e3 * ones(3, 1); M = diag(m); k = 2e6 * ones(3, 1); K = zeros(3, 3); K(1, 1) = k(1) + k(2); K(1, 2) = -k(2); K(2, 1) = -k(2); K(2, 2) = k(2) + k(3); K(2, 3) = -k(3); K(3, 2) = -k(3); K(3, 3) = k(3);

阻尼采用Rayleigh阻尼，即C = α·M + β·K。这是一种工程上常用的简化方式，好处是保证阻尼矩阵的正定性，而且在频域里有明确的物理含义。取α=0.3，β=2×10^(-4)，这样三阶模态的阻尼比大约在1%到2%之间，比较接近真实钢结构或混凝土结构的水平。

有了M、C、K，状态空间模型可以按标准的"位移-速度"状态向量写出：

n = 3; A_ss = [zeros(n), eye(n); -M\K, -M\C]; B_ss = [zeros(n); inv(M)]; C_ss = eye(6); % 观测全部位移和速度 D_ss = zeros(6, n); sys = ss(A_ss, B_ss, C_ss, D_ss);

状态向量是6维的，前3维是各层位移，后3维是各层速度。C_ss取单位阵意味着我们理论上可以观测所有状态量，但在实际工程中通常只测位移或加速度响应，这里为了模拟真实情况，只取位移输出（前3列）用作后续识别。

先算一下理论上的模态参数作对照基准，这很重要，只有知道"正确答案"才能验证识别结果的准确性：

[V, D] = eig(K, M); % 广义特征值问题 omega2 = diag(D); fn_theory = sqrt(omega2) / (2 * pi); % 排序后即为结构前三阶固有频率

把模型参数代进去能算出前三阶固有频率分别约为1.45 Hz、4.06 Hz、5.93 Hz。

3.2 环境激励下的响应模拟

环境激励的本质是宽频随机激励，工程上常用高斯白噪声近似。下面用lsim对系统施加白噪声激励并求响应：

fs = 200; % 采样频率 200 Hz，覆盖前两阶频率绰绰有余 dt = 1 / fs; T_dur = 60; % 仿真时长 60 秒 t = (0:dt:T_dur-dt)'; rng(1); % 固定随机种子，保证结果可复现 u = 1e4 * randn(length(t), 3); % 三个质量块各自受白噪声激励 u(:, 2:3) = 0.5 * u(:, 2:3); % 主激励放在第一层，模拟基底激励 [y, ~] = lsim(sys, u, t); y_disp = y(:, 1:3); % 只取前3列位移响应

这里激励幅值取1e4是为了得到幅值在10^(-3)量级的位移响应，更接近实际结构的微振水平。采样频率选200Hz是因为最高关注频率不到6Hz，200Hz已经留足余量，能看清高频区噪声的分布情况。

实际测量信号永远伴随噪声，为了检验SSI-COV的抗噪能力，给响应叠加上信噪比30dB的高斯白噪声：

y_noisy = awgn(y_disp, 30, 'measured');

到这里，仿真数据就准备好了。

4. Matlab实现SSI-COV全流程

4.1 主函数设计

我习惯把SSI-COV的核心流程封装成一个函数，输入是响应数据、采样频率、Hankel块数和系统阶数，输出是识别出的频率、阻尼比和振型。函数主体如下（这是完整可直接运行的版本）：

function [fn, zeta, phi, A_rec] = ssi_cov(y, fs, i_blocks, order) % SSI-COV 协方差驱动随机子空间识别 % 输入： % y - 响应数据，nCh x N，行为测点通道，列为时间采样 % fs - 采样频率，Hz % i_blocks - Hankel矩阵的过去/未来块数 i % order - 系统阶数，通常为 2*模态数 % 输出： % fn - 识别频率，单位 Hz % zeta - 阻尼比，无量纲 % phi - 复振型矩阵，每一列对应一阶模态 % A_rec - 识别出的状态矩阵 [nCh, N] = size(y); % 1. 构造分块Hankel矩阵 H = zeros(2 * i_blocks * nCh, N - 2 * i_blocks + 1); for k = 1:2 * i_blocks H((k-1)*nCh+1 : k*nCh, :) = y(:, k : k+N-2*i_blocks); end Yp = H(1 : i_blocks*nCh, :); % 过去块 Yf = H(i_blocks*nCh+1 : end, :); % 未来块 % 2. 构造Toeplitz协方差矩阵并标准化 T = (Yf * Yp') / size(Yp, 2); % 3. 对Toeplitz矩阵做SVD分解 [U, S, V] = svd(T, 'econ'); if order > size(S, 1) error('order 超过SVD分解得到的秩，请减小order或增大Hankel块数'); end U1 = U(:, 1:order); S1 = S(1:order, 1:order); % 4. 还原可观测性矩阵和输出矩阵C Oi = U1 * sqrt(S1); C_rec = Oi(1:nCh, :); % 5. 利用位移结构最小二乘求解系统矩阵A A_rec = pinv(Oi(1:(i_blocks-1)*nCh, :)) * Oi(nCh+1:i_blocks*nCh, :); % 6. 特征值分解并转换到连续时间域 [Psi, Lambda] = eig(A_rec); lambda_d = diag(Lambda); % 离散时间特征值 s_c = log(lambda_d) / dt_act; % 连续时间特征值 dt_act = 1 / fs; % 7. 提取频率、阻尼比、振型 fn = abs(s_c) / (2 * pi); zeta = -real(s_c) ./ abs(s_c); phi = C_rec * Psi; % 振型矩阵 end

细心的读者会发现，第3步SVD截断依赖order，而实际工程中我们并不知道order该取多少。这正是稳定图要解决的问题，后面我会单独展开。上面的函数先给固定order时用的版本。

4.2 从Toeplitz矩阵到SVD截断的细节

这里有三处容易被忽略但决定成败的细节：

第一，Hankel矩阵的列数。我在代码里用的是N-2i_blocks+1，这是为了保证过去块和未来块都有完整数据。有些实现会直接用N-2i_blocks（少一列），误差很小，但要在代码里保持统一。

第二，Toeplitz矩阵计算时的归一化。严格推导时，T的元素是互协方差E[y(k+i)·y(k)^T]，所以除以采样点数。工程上除以N或N-1差别不大，但除以多少必须和后面的SVD结果一起理解，因为奇异值幅值会随缩放变化。

第三，SVD截断后系统矩阵A的维度是order×order。A的特征值如果出现实部为正的极点，在物理上是模态失稳的现象，说明识别出了虚假模态（真实结构不可能发散），在筛选时可以剔除。

4.3 模态参数提取：离散特征值到连续频率阻尼

这一步是SSI-COV误差最容易被放大的地方。随机子空间识别得到的A是离散时间状态矩阵，它的特征值是离散时间极点的形式，需要转换到连续时间极点。

离散特征值λ_d与连续特征值s的关系是λ_d = e^(s·Δt)，反过来：

s = ln(λ_d) / Δt

s通常是复数，实部对应衰减率，虚部对应有阻尼振动频率。系统的无阻尼固有频率、有阻尼固有频率和阻尼比之间有如下数学关系：

s = -ζ·ω_n ± j·ω_d |s| = ω_n ω_d = Im(s) ζ = -Re(s) / |s|

所以：

fn = abs(s_c) / (2 * pi); % 无阻尼固有频率 fd = imag(s_c) / (2 * pi); % 有阻尼频率 zeta = -real(s_c) ./ abs(s_c);

工程上阻尼比小于10%的时候，f_n和f_d相差不到0.5%，报告频率用哪个都行，但要注明。阻尼比直接由实部和模的比值算出，不需要额外的数值微分，这也是SSI-COV测阻尼比比较准的原因之一。

振型提取稍特殊一点：振型矩阵的每一列是特征方程中ψ的列，但测点只关注输出位置。由于输出向量y(k)=C·x(k)，对应第r阶模态的输出振型为C·ψ_r，也就是phi = C_rec * Psi的每一列。

4.4 一次完整识别：结果长什么样

拿前面生成的仿真数据去做固定阶数识别。取i_blocks=10，order=6（3阶模态），运行主流程得到：

i_blocks = 10; order = 6; [fn_ssi, zeta_ssi, phi_ssi] = ssi_cov(y_noisy, fs, i_blocks, order);

稳定后识别的结果大致如下（由于每次白噪声序列不同会有轻微浮动，但趋势一致）：

频率误差普遍在0.5%以内，阻尼比误差在10%以内，对信噪比30dB的数据来说，这个精度已经相当能打。我拿这个方法处理过实测的桥梁微振数据和楼板振动数据，频率识别的精度足够支撑后续模型修正工作；阻尼比虽然比频率容易飘，但在现场条件下比频域法给出的结果稳定得多。

5. 稳定性图：让系统阶次自己"说话"

5.1 为什么要引入稳定图

SSI-COV在工程应用中最实际的问题就是定阶——到底保留多少个系统阶次合适？奇异值曲线能给出参考，但实际数据信噪比不足时，叠加虚假模态的情况几乎难以避免。这让固定order的识别方法显得有点脆弱：取大了引入大量虚假极点和噪声极点，取小了又可能漏掉弱激励的真实模态。

稳定图的思路很朴素：不管真实模态还是虚假模态，都会随着系统阶次变化而移动。但真实模态在所有阶次下都会保持稳定，而虚假模态则忽隐忽现、参数漂来漂去。把不同阶次下的识别极点的频率、阻尼、振型画在同一张图上，保留那些"稳如泰山"的极点，就是物理模态。

5.2 稳定图的三个判据

工程上最常用的稳定判据有三个：

• 频率稳定性：相邻阶次下频率偏差小于1%，即|f_i - f_j| / f_i < 0.01
• 阻尼稳定性：相邻阶次下阻尼比差值的绝对值小于一定阈值，常用0.05（绝对值）或10%（相对值）
• 振型稳定性：两个振型之间的MAC值大于0.95

MAC（Modal Assurance Criterion，模态置信准则）是衡量两个振型关联度的经典指标：

MAC = |ψ_i^H · ψ_j|² / (|ψ_i^H · ψ_i| · |ψ_j^H · ψ_j|)

MAC越接近1说明两个振型越一致，大于0.95就认为振型稳定。

5.3 稳定图绘制的Matlab实现

画稳定图需要循环不同order，把每次识别出的频率、阻尼、MAC整理成表，再标定稳定点：

% 准备存储 orders = 2:2:40; % 扫描系统阶次从2到40，步长2 n_orders = length(orders); ptr = 1; % 每个阶次做一次SSI-COV for k = 1:n_orders order = orders(k); [fn_o, zeta_o, phi_o] = ssi_cov(y_noisy, fs, 10, order); % 只保留实部为负的物理极点（剔除发散极点） keep = real(s_c_here) < 0; % 需要在ssi_cov中额外返回s_c，或重写函数 fn_o = fn_o(keep); zeta_o = zeta_o(keep); phi_o = phi_o(:, keep); % 记录到单元格中 freq_all{k} = fn_o; zeta_all{k} = zeta_o; phi_all{k} = phi_o; end

实际上要完整实现稳定图，我建议对ssi_cov函数返回连续特征值s_c和振型phi，后续再统一筛选。画图的核心是把每个阶次算出的频率按对应纵坐标order值画一个点，稳定点用不同颜色标记：

figure; hold on; for k = 1:n_orders y_coord = orders(k) * ones(size(freq_all{k})); % 先画所有点 plot(freq_all{k}, y_coord, 'k.', 'MarkerSize', 4); end xlabel('频率 (Hz)'); ylabel('系统阶次 (2n)'); ylim([0 max(orders)+2]);

叠加稳定点判断逻辑后，图上会出现若干条竖直的"柱子"——柱子密集的位置就是真实模态的频率位置。

5.4 实际操作中怎么读稳定图

我扫完order 2到40后，稳定图上有三处竖线最明显，分别在1.45 Hz、4.06 Hz和5.93 Hz附近，正好对应三阶真实模态。有些高order区域偶尔冒出其他频率的散点，但它们在相邻阶次之间到处飘，或者MAC值达不到0.95，直接忽略。

一个我自己用了很久的习惯：先把频率稳定和极点发散条件作为硬性筛选，阻尼比作为软性参考。阻尼比判据容易被噪声干扰，我会把阻尼阈限放宽一点，否则很多真实模态因为阻尼不稳定被误杀。频率和振型才是判断模态真实性的可靠指标。

6. 实操过程中的坑与排查方法

6.1 数据预处理：去均值和趋势项

SSI-COV本质上是基于协方差统计的方法，如果数据里存在非零均值或者缓慢漂移的趋势项，协方差矩阵会被低频成分污染，导致识别出接近0 Hz的虚假极点。实测数据尤其容易遇到这个问题，传感器温漂、线缆干扰会造成基线漂移。

我每次拿到数据的第一件事是减均值，这是必须的。如果信号有明显趋势项，再做一个高通滤波或多项式去趋势：

y_disp = detrend(y_disp, 'linear'); % 去除线性趋势

有些现场信号还有明显工频干扰（50Hz或60Hz），SSI-COV会把工频成分当作一个极稳定的"伪模态"识别出来。工频确实稳定，稳定图上会形成非常干净的竖线，但它不是结构模态。所以信号通过低通滤波器、把截止频率设在关注频段以上的1.2倍，是必要的预处理步骤。

6.2 Hankel块数i_blocks怎么选

i_blocks这个参数直接影响Toeplitz矩阵的大小和数据的利用效率。选小了，协方差信息不足，低阶模态可能识别不出来；选大了，Toeplitz矩阵变得很庞大，计算慢了不说，还容易把更多噪声细节纳入模型中，带来大量虚假极点。

我的一般经验是三个原则：

• i_blocks2采样时长要覆盖至少20个最低关注周期的数据
• i_blocks取数据总采样点数的5%到10%之间，通常10到30之间足够
• 在数据充裕的情况下尽量多试几组取值，看结果是否稳定

对上面的仿真数据，总采样点数12000，i_blocks取10时频率结果已经很稳定，取20时结果几乎不变，取5时频率略有偏差。可见这个参数有较宽的合理区间，不太需要严格寻优。

6.3 为什么识别出的频率总是偏低一点点

用SSI-COV识别出的阻尼比偏正、频率偏负，是很多初学者会注意到的细节。这不是代码bug，而是离散化采样的固有偏差：时域识别方法本质上是把连续系统映射到离散时间域，采样率不够高时，这种映射存在系统偏差。要减小偏差，就得提高采样频率Δt。一般保证最高关注频率对应的每周期采样点数不小于10个，偏差可以忽略。

频率偏差还有一个来源是数据长度不够。SSI-COV利用Toeplitz矩阵估计协方差，数据越短，协方差估计的方差越大，识别结果越容易系统性偏小。60秒数据对1.5Hz的结构意味着约90个振动周期，这是比较稳妥的底限；如果现场只能采到20秒数据，就需要接受精度损失。

6.4 虚假模态的识别与剔除

虚假模态是SSI-COV绕不开的话题。除稳定图外，我总结出几个很实用的筛选经验：

• 剔除极点实部为正的发散伪模态
• 剔除阻尼比小于0%或大于15%的极点，真实结构的阻尼比极少超过10%
• 剔除阻尼比恰好落在边界上的可疑极点，真实模态的阻尼比一般不会奇怪地取整
• 振型MAC值小于0.7的极点不要信，欠激励的模态振型相关性很差

大量的实践中，"三连筛"——频率稳定、阻尼在合理范围、振型MAC达标——能过滤掉九成以上的虚假模态。剩下零星几个漏网的，结合振型形状是否满足结构物理规律（比如梁的高阶振型零点位置是否合理）做人工判断。

7. 一点个人体会

SSI-COV这套方法我前前后后用了快三年，从一开始只会在仿真数据里跑通，到后来处理实测信号，最大的感受是：方法本身很优雅，但工程落地全在细节。仿真数据里随便怎么调参都能识别得漂漂亮亮，实测数据一到手，去均值、去趋势、滤波、选块数、判稳定，每一步都可能让结果翻车。现阶段如果你还需要一个（起点足够好的）实现方案，上面这套代码和流程能直接跑通三自由度到几十自由度的线性结构。这个内容后续还可以扩展的方向也很多，比如结合频域分解法做交叉验证、扩展成SSI-DATA对比计算效率、或者引入自动聚类算法让稳定图全自动判读，都是价值很高的延伸方向。