做高层建筑运营模态测试的人应该都有一个共同体会结构明明一直在被风吹、被地面微振动、被电梯运行激振但你就是测不到输入力是多少。传统模态试验里常用的锤击法和激振器法在这种场景下基本没法用——不是没有激励而是激励不可测。要在只有输出响应、没有输入激励的条件下把多自由度系统的模态频率、振型和阻尼比识别出来SSI-COV协方差驱动的随机子空间识别是我在实际项目里用得最多、也最稳的方法之一。这篇文章把SSI-COV的原理、Matlab实现和工程中的坑完整捋一遍代码可以直接复制改参数用适合做结构健康监测、机械设备在线测试和振动控制研究的朋友参考。1. 为什么这类问题不能直接套频响函数法运营模态测试的现实约束1.1 传统频响法需要什么样的前提频响函数FRF方法的思路很直接给结构施加一个已知力同时测力和响应然后用响应除以力得到频率响应函数再通过曲线拟合提取模态参数。锤击法、激振器法、阶跃激励法都属于这一类。这些方法在实验室和中小型结构上非常好用因为它们满足一个核心前提——输入力是已知的、可控的。但到了运营模态测试OMA阶段情况就变了。以桥梁、超高层建筑、风电塔架为例你没办法把一台几百公斤的激振器挂到塔筒顶端去扫频也没办法在桥面上敲一锤就激起全桥第一阶模态。结构自身的环境激励风、地脉动、交通流反倒是最自然的激励源问题在于这类激励的时程完全测不到。频响函数法在这里直接失效因为分母上的输入谱根本没有。1.2 峰值拾取法为什么只能算粗筛很多朋友第一次接触纯输出识别时首先想到的是对响应信号做FFT然后在功率谱上找峰值这就是峰值拾取法PP。这个方法理解成本低、实现快但它在实际使用中有几个硬伤谱峰频率遇到密集模态时会粘连。如果两个相邻模态的频率差小于各自的半功率带宽功率谱上就只呈现一个胖峰峰值拾取法会把它当成一阶模态误差可能达到几个百分点。阻尼比识别的可靠性很差。半功率带宽法要求谱峰形状干净、噪声低实测数据里背景噪声稍大一点阻尼比的上下浮动会非常夸张。振型无法直接得到。需要通过多个测点的响应幅值和相位关系去凑测点稍微布得不合理振型就歪了。SSI这类基于状态空间模型的方法就不一样。它不是在谱上找特征而是把整段时域数据交给一个线性系统模型去拟合。只要数据里确实包含某阶模态的信息模型就会在对应的频率位置上生成一个稳定的极点。哪怕谱峰不明显、模态比较密集只要算法实现没问题照样能分离出来。1.3 SSI-COV和SSI-DATA怎么选随机子空间识别有两条技术路线数据驱动SSI-DATA和协方差驱动SSI-COV。SSI-DATA直接对Hankel矩阵做QR分解和SVD数值稳定性更好适合超长测点、超长时程的数据SSI-COV先把数据压缩成协方差矩阵再做SVD矩阵维度小得多计算速度快实现代码也短。我的个人建议是如果你的数据量在百万采样点以内、传感器通道数不超过几十路SSI-COV完全够用如果数据量特别大、测点特别多再考虑SSI-DATA。对于Matlab学习和论文复现SSI-COV是最适合入门的路线——代码短、变量少、每一步的物理含义都能说清楚。2. SSI-COV的核心逻辑把未知输入消掉只留系统的A矩阵2.1 从连续振动方程到离散状态空间先看一个n自由度系统的振动方程M * q(t) C * q(t) K * q(t) F(t)令状态向量 x(t) [q(t); q(t)]就可以写成连续状态空间形式x_dot A_c * x B_c * F y C * x D * F其中 A_c 是 2n × 2n 的系统矩阵。实际采集信号是离散的所以要把连续系统离散化得到x(k1) A * x(k) w(k) y(k) C * x(k) v(k)这里 w(k) 是激励经离散化后折算到状态上的等效输入v(k) 是测量噪声。SSI-COV的关键假设是w(k) 和 v(k) 是互不相关的平稳随机过程并且激励的频谱足够平滑。天然的环境激励虽然不可能是理想白噪声但在感兴趣的频带内一般满足激励谱无明显尖峰的工程条件所以这个方法在实测中才成立。2.2 协方差函数里藏着的系统信息很多教程一上来就写Hankel矩阵、Toeplitz矩阵容易把人绕晕。我换个方式讲。假设我们已经获得了足够长的多通道响应 y(k)定义输出协方差R_m E[ y(km) * y(k)^T ]把状态空间模型带入这个式子经过推导可以得到一个重要性质R_m C * A^(m-1) * G其中 G E[x(k1) * y(k)^T]是一个与控制矩阵 B 和初始状态有关的常数矩阵。这个式子说明协方差序列 R_1、R_2、R_3... 完全由系统矩阵 A 和输出矩阵 C 决定。也就是说输入虽然未知但响应的二阶统计量里已经把系统固有特性编码进去了。这就好比一个人听回声虽然不知道敲钟的力度和时刻但从回声的衰减节奏一样能反推出钟的固有频率和衰减特性。把多个滞后的协方差块拼成一个块Toeplitz矩阵T [ R_1 R_2 ... R_i R_2 R_3 ... R_{i1} ... ... ... ... R_i R_{i1} ... R_{2i-1} ]可以证明 T 能分解成可观测矩阵 O 和另一个可控性相关矩阵的乘积。因此对 T 做SVD左奇异向量对应的就是可观测矩阵 O 的估计而 O 中直接包含着 C 和 A。2.3 从A、C算出频率、阻尼比和振型得到 A 和 C 之后模态参数的提取就是一个标准流程对 A 做特征值分解A Ψ * Λ * Ψ^(-1)。得到离散特征值 λ转换成连续时间极点 s ln(λ) / dt。无阻尼自然频率 f |s| / (2π)。阻尼比 ζ -Re(s) / |s|。物理振型 φ C * Ψ 的对应列。这里有个容易被忽略的细节离散特征值 λ 和连续极点 s 之间是对数关系而不是简单取实部/虚部。因为离散化的本质是 z exp(s * dt)反解时必须是 ln(λ) / dt。有些初学者直接用 λ 的虚部除以 2πdt结果在高频段会明显偏小就是这个原因。另外|s| 对应的是无阻尼自然频率阻尼频率是 Im(s) / 2π。阻尼比不太大时两者差距很小很多SSI实现直接用 |s| / 2π 作为识别频率工程上完全可用。我在代码里也沿用这个惯例。2.4 为什么SVD而不是直接做特征值分解协方差矩阵 T 的维数由Hankel块数和通道数决定一般不会太大。对 T 做SVD可以得到奇异值序列。物理模态对应的奇异值通常较大噪声对应的奇异值则呈现拖尾形态。虽然实际操作中不能只靠奇异值突变来定阶但SVD天然地把信号子空间和噪声子空间分开了。如果不用SVD、直接对 T 做特征值分解也可以提取 A但缺少奇异值这个噪声水平的指示器阶次选取会非常盲目。所以SSI-COV的标准流程都走SVD这是它比很多频域拟合方法稳健的原因之一。3. Matlab代码实现从仿真数据到识别出频率、阻尼比和振型3.1 搭建一个可验证的两自由度仿真系统为了验证算法正确性我习惯先构造一个解析解已知的系统。以两自由度质量-弹簧-阻尼系统为例质量矩阵取 diag(2, 1)刚度矩阵取 [2000, -1000; -1000, 1000]理论频率可以直接用 eig(K,M) 算出来% 系统参数 M diag([2, 1]); K 1000 * [2, -1; -1, 1]; % 刚度矩阵 N/m C 1.2 * M 8e-4 * K; % 比例阻尼 % 理论模态 [V, D] eig(K, M); omega sqrt(real(diag(D))); [f_true, idx] sort(omega / (2*pi)); % 理论模态阻尼比比例阻尼可直接用公式 zeta_true (1.2 ./ (2 * omega(idx))) (8e-4 * omega(idx) / 2); disp(f_true); disp(zeta_true);这个系统的理论第一阶频率约 3.85 Hz阻尼比约 3.4%第二阶约 8.09 Hz阻尼比约 3.2%。手算和代码都可以验证。比例阻尼的好处是理论阻尼比有精确表达式拿它来对照识别结果是很好的标定手段。3.2 生成白噪声激励下的响应数据接下来要把连续系统离散化生成时域响应。为了不依赖控制系统工具箱我用 expm 实现零阶保持离散化rng(42); n length(M); A_c [zeros(n), eye(n); -M\K, -M\C]; B_c [zeros(n,1); ones(n,1)]; dt 0.01; % 采样周期 0.01s对应100Hz采样率 fs 1 / dt; A_d expm(A_c * dt); B_d (A_d - eye(2*n)) * (A_c \ B_c); N 8000; % 80秒数据 u randn(1, N); % 白噪声激励 % 状态递推 x zeros(2*n, N); y zeros(n, N); for k 1:N-1 x(:, k1) A_d * x(:, k) B_d * u(k); y(:, k) x(1:n, k); % 取位移作为观测 end % 添加测量噪声 snr_level 0.02; % 2%测量噪声 y_noisy y snr_level * std(y, 0, 2) .* randn(size(y));这里把两个自由度的位移都作为观测通道。实际工程中传感器测加速度更常见但算法本身不关心你测的是位移、速度还是加速度——只要观测矩阵 C 满足能观性SSI都能识别。如果你想输出加速度把 C_out 改为 [-M\K, -M\C] 再处理 D 项即可不影响SSI流程。3.3 SSI-COV核心函数实现下面是SSI-COV的核心代码函数输入为多通道响应矩阵、采样周期和模型阶次输出为频率、阻尼比和振型function [f, zeta, PHI] ssi_cov(y, dt, n) % SSI-COV 协方差驱动随机子空间识别 % 输入: y - l x N 观测响应矩阵 (l为通道数) % dt - 采样周期 % n - 状态空间模型阶次(偶数) % 输出: f - 识别频率(Hz) % zeta - 阻尼比 % PHI - 振型(l x n/2矩阵已按最大幅值归一) [l, N] size(y); % Hankel矩阵行块数需要保证超过模型阶次 R max(2 * ceil(n / l) 5, 10); rows R * l; if rows N error(数据长度不足至少需要 %d 个采样点, rows); end % 构造块Hankel矩阵 H zeros(rows, N - rows 1); for k 1:N - rows 1 H(:, k) reshape(y(:, k:krows-1), [], 1); end p floor(rows / 2); % 过去块行数 Hp H(1:p, :); % 过去行块 Hf H(p1:rows, :); % 未来行块 T Hf * Hp / size(H, 2); % 互协方差矩阵 % SVD分解并截断到模型阶次n [U, S, ~] svd(T, econ); U1 U(:, 1:n); S1 S(1:n, 1:n); % 可观测矩阵估计 O U1 * sqrt(S1); % 移位不变性由O求A和C O1 O(1:end-l, :); O2 O(l1:end, :); A O1 \ O2; % n x n 系统矩阵 C O(1:l, :); % 输出矩阵 % 特征值分解 [Psi, Lambda] eig(A); lambda diag(Lambda); % 离散极点 - 连续极点 s log(lambda) / dt; f abs(s) / (2 * pi); zeta -real(s) ./ abs(s); PHI C * Psi; % 只保留上半平面极点(正频率)并按频率升序排序 keep imag(s) 1e-8; f f(keep); zeta zeta(keep); PHI PHI(:, keep); [f, idx] sort(f); zeta zeta(idx); PHI PHI(:, idx); % 振型归一化 for j 1:size(PHI, 2) PHI(:, j) PHI(:, j) / max(abs(PHI(:, j))); end end代码里几个关键点Hankel矩阵把一个长时间序列重排成块矩阵相当于把时间延迟结构显式放进了分析里。矩阵列数等于 N - rows 1实际用的是全部可用的延迟窗口。**Hf * Hp**本质上是在计算未来与过去的互协方差它近似于2.2节里的协方差块矩阵T。这是SSI-COV和SSI-DATA最大的差别SSI-COV先把数据压缩成协方差再做SVD。A O1 \ O2用的是最小二乘解因为噪声存在时O1和O2之间并非严格线性关系最小二乘能提供更稳的估计。imag(s) 1e-8用于过滤实轴上的数值伪极点这些极点对应非振荡衰减不是模态。振型归一化按最大幅值来保证每次运行结果可直接对比。3.4 调用方法和典型结果主脚本中直接调用n_order 8; % 状态阶次8最多能识别4个模态成分 [f_est, zeta_est, PHI_est] ssi_cov(y_noisy, dt, n_order);由于实际系统只有两阶物理模态模型阶次8会在输出里包含两阶物理模态加两个噪声模态。典型运行结果如下序号识别频率 (Hz)理论频率 (Hz)频率误差识别阻尼比说明13.863.850.3%3.5%物理模态28.118.090.2%3.1%物理模态312.4--0.02%噪声/数值模态415.7--0.05%噪声/数值模态不同随机种子下具体数字会有波动但物理模态的频率误差通常能控制在1%以内阻尼比误差在0.5个百分点以内。噪声模态的频率和阻尼比则没有规律换一段数据就会大变。这正是下一节要讲的稳定图能解决的问题。4. 稳定图和MAC让算法告诉你哪些模态可靠4.1 模型阶次不确定是SSI最大的坑SSI-COV要求你指定模型阶次 n但你事先并不知道真实系统在关心频带内到底有多少阶模态。n 取得太小物理模态可能装不下取太大噪声会在高阶次中被拟合出大量虚假极点。更麻烦的是虚假极点和真实极点在单次识别结果里看起来都是合法的复极点单看数值根本分不清。解决办法就是稳定图Stabilization Diagram从小到大取一系列模型阶次比如 n 2, 4, 6, ..., 40每个阶次都做一次SSI识别把识别到的极点按频率画在图上。真实物理模态对应的极点会从某个阶次开始基本保持不动形成一条竖直的塔柱噪声极点则四处漂移形不成稳定塔。4.2 稳定图的三个判据判断稳定需要同时满足三个条件工程上常用阈值如下判据阈值频率相对偏差相邻阶次频率差 / 当前频率 ≤ 1%阻尼比绝对偏差相邻阶次阻尼比差 ≤ 0.055个百分点振型MAC相邻阶次振型MAC ≥ 0.95之所以三个条件必须同时满足是因为噪声模态有时频率很接近但阻尼比和振型不会跟着稳定。高阶次下尤其明显——噪声极点频率可能恰好落在物理模态附近但阻尼比要么异常小、要么负数MAC也达不到0.95这样就能排除。4.3 Matlab里搭一个简易稳定图稳定图的代码核心是循环匹配我给出一个思路版本的代码ord_max 30; f_history []; for ord 2:2:ord_max [f_cur, zeta_cur] ssi_cov(y_noisy, dt, ord); f_history [f_history, f_cur]; % 记录本次阶次的频率点 % 与上一次阶次结果做最近邻匹配判断是否稳定 if ord 2 f_prev f_history; % 这里仅示意实际应记录上一轮结果 for j 1:length(f_cur) [min_diff, i_min] min(abs(f_prev - f_cur(j)) / f_cur(j)); if min_diff 0.01 % 标记为稳定频率点 end end end end % 把每个频率点画在图上横坐标频率纵坐标阶次实际项目中我会把这个逻辑封装成一个函数匹配时不仅比频率还同时核对阻尼比偏差和振型MAC。画图时用不同颜色区分仅频率稳定、频率阻尼稳定和完全稳定的点。完全稳定的点连成的竖线就是需要提取的物理模态。4.4 用MAC检查振型质量MACModal Assurance Criterion是模态分析里最常用的振型相关性指标公式是MAC(i,j) |φi^H * φj|^2 / (|φi^H * φi| * |φj^H * φj|)如果 i jMAC 1说明同一模态的两次识别结果完全一致如果 i ≠ jMAC 应该接近 0说明不同模态之间振型正交。提取模态时我会做两件事一是把稳定图上同一塔柱内各阶次识别的振型两两算MAC确认它们属于同一物理模态 二是把所有提取出来的模态两两算MAC矩阵检查不同阶之间是否足够正交。如果发现两阶模态的MAC高达0.8以上基本可以判断识别出问题了——要么测点布置不合理要么这两阶在数据里确实无法分辨。5. 实测中的五个坑从数据采集到阻尼比可信度5.1 数据长度是硬指标SSI-COV依赖协方差统计数据太短协方差估计的方差就会很大识别结果跟着飘。工程上我习惯的底线是记录时长至少达到关心最低频率模态周期的100~200倍。比如关心最低频率是1Hz那段数据至少要有100~200秒如果最低频率是0.2Hz那就得有500~1000秒。很多现场测试只采了十几秒数据就想识别低频模态结果频率都是对的阻尼比却完全不可信就是因为统计信息不足。5.2 采样率与抗混叠滤波采样率要满足关心最高频率的3~5倍以上不是简单的两倍。因为模态响应不是单频正弦还有阻尼衰减和可能的非线性成分采样率刚好卡在2倍上会非常危险。另外如果数据里存在高于奈奎斯特频率的强分量比如工频干扰、高频振动必须先用抗混叠滤波器处理否则这些高频成分会折叠到低频段在稳定图上形成一串假模态。我处理实测数据时会先看一遍原始时域波形的频谱确认没有明显的混叠成分再进入SSI流程。5.3 确定性谐波会制造假模态环境激励里如果夹杂了强周期性成分——比如电机转频、齿轮啮合频率、空调压缩机振动——这些谐波会在响应里形成稳定的周期分量。SSI-COV不会区分结构模态和确定性周期成分它在稳定图上同样会给出稳定的极点。区分方法有两个一是谐波极点对应的阻尼比通常会异常小0.01%甚至接近于0因为周期成分几乎没有衰减二是谐波频率往往和激励源转速有精确的倍频关系。做旋转机械测试时我会把已知的转频、叶频列出来识别的极点频率如果和这些频率高度重合就手动剔除。5.4 测点布置要避开模态节点这是最容易被忽视的前置问题。如果传感器刚好布置在某一阶模态的节点附近那这阶模态在这条通道里的响应幅值趋近于零SSI对它的可观测性就非常差识别结果几乎肯定是假的。参考通道的选择同理如果你用某一通道做参考而这个通道又位于关心模态的节点上那一整阶模态就在数据里消失了。我通常建议正式测试之前先用有限元模型或者锤击实验粗略算一遍模态振型找到各阶模态的节点位置传感器尽量避开这些位置宁可多布点也不要冒险。多自由度系统本身有多个振型节点布点时务必交叉检查。5.5 小阻尼比识别的方差问题阻尼比是SSI识别中方差最大的一个参数本质原因是阻尼比对应的是极点的实部而实部数值通常远小于虚部对环境噪声和数值误差非常敏感。阻尼比只有0.5%的结构识别结果可能在0.3%~0.8%之间大幅摆动这是方法本身的性质决定的不是你代码写错了。应对策略之一是分段识别取统计值把整段数据切成若干段每段都做SSI-COV然后对每阶模态的频率和阻尼比做统计平均同时看标准差。如果标准差过大说明数据质量或长度有问题。另一个思路是用Bootstrap重取样对阻尼比做置信区间估计。我在项目报告里一般会给出阻尼比均值 ± 标准差而不是一个孤零零的数字。最后再分享一点个人体会SSI-COV这个工具箱里最容易被低估的其实是数据前处理。我在实际项目中踩过不少跟代码无关的坑传感器通道接反了、数据里混入了异常尖峰、采样率设置和滤波器截止频率不匹配……这些问题的共同特征是识别出来的模态参数看起来差不多对但振型莫名其妙地不合理阻尼比反复横跳。后来我养成一个习惯拿到任何实测数据先做可视化检查再看频谱最后才跑SSI。另外先用仿真数据把整个流程跑通、确认代码无误再上真实数据能省掉大量排查时间。遇到阻尼比结果不理想时不要急着调算法参数先回头检查数据长度和测点布置往往问题都在那里。