我前阵子帮一个做土木的朋友处理过一串实测楼板振动数据慢慢意识到一个问题很多搞结构的人手里晃着好数据却对怎么把模态参数干净地摘出来这一步犯怵。频域峰值法简单但主观拟合传递函数又对激励有硬性要求到了现场实测这种环境激励占主导的场景传统方法常常翻车。这次就把我在多自由度系统上用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(k1) 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下半部分叫未来YfYp [y(0) y(1) ... y(N-2i)] [y(1) y(2) ... y(N-2i1)] ... [y(i-1) y(i) ... y(N-i-1)]Yf [y(i) y(i1) ... y(N-i)] [y(i1) y(i2)... y(N-i1)] ... [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 · Γ_iO_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(L1:iL, :)第五步对A做特征值分解从特征值里解出频率和阻尼比。这一套流程逻辑非常严密每一步都有扎实的线性代数支撑这也是它比频域峰值法更让人放心的原因。2.3 几个关键点的直观理解这里用大白话解释一下为什么SVD能帮我们定阶。想象T矩阵的信息分为真实振动和噪声两部分前者在奇异值谱上表现为较大的奇异值后者散布在小奇异值上。做SVD相当于把矩阵按信息重要程度重新分解排序只要真实模态的奇异值明显大于噪声奇异值截断点就一目了然。实际数据如果信噪比不高奇异值可能平滑衰减没有明显台阶这时候就需要稳定图辅助判断后面专门讲。另一个关键点是每阶物理模态对应一对共轭复数极点所以在状态空间里系统阶次是2n而不是n。代码里如果扫描到5阶模态对应系统阶数应该设为10。新手经常在这里犯迷糊。3. 多自由度系统仿真数据准备3.1 三自由度剪切结构模型的建立为了验证SSI-COV的实现效果我建了一个经典的三层剪切型结构模型模型简图就是三个质量块串联层间刚度和阻尼集中。取质量m1m2m35000 kg层间刚度k1k2k32×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是因为最高关注频率不到6Hz200Hz已经留足余量能看清高频区噪声的分布情况。实际测量信号永远伴随噪声为了检验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)*nCh1 : k*nCh, :) y(:, k : kN-2*i_blocks); end Yp H(1 : i_blocks*nCh, :); % 过去块 Yf H(i_blocks*nCh1 : 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(nCh1: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_blocks1这是为了保证过去块和未来块都有完整数据。有些实现会直接用N-2i_blocks少一列误差很小但要在代码里保持统一。第二Toeplitz矩阵计算时的归一化。严格推导时T的元素是互协方差E[y(ki)·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) / Δts通常是复数实部对应衰减率虚部对应有阻尼振动频率。系统的无阻尼固有频率、有阻尼固有频率和阻尼比之间有如下数学关系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_blocks10order63阶模态运行主流程得到i_blocks 10; order 6; [fn_ssi, zeta_ssi, phi_ssi] ssi_cov(y_noisy, fs, i_blocks, order);稳定后识别的结果大致如下由于每次白噪声序列不同会有轻微浮动但趋势一致模态阶数理论频率(Hz)SSI-COV识别频率(Hz)理论阻尼比(%)SSI-COV识别阻尼比(%)11.4521.4511.271.3124.0684.0660.961.0235.9365.9310.840.91频率误差普遍在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.95MACModal 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或60HzSSI-COV会把工频成分当作一个极稳定的伪模态识别出来。工频确实稳定稳定图上会形成非常干净的竖线但它不是结构模态。所以信号通过低通滤波器、把截止频率设在关注频段以上的1.2倍是必要的预处理步骤。6.2 Hankel块数i_blocks怎么选i_blocks这个参数直接影响Toeplitz矩阵的大小和数据的利用效率。选小了协方差信息不足低阶模态可能识别不出来选大了Toeplitz矩阵变得很庞大计算慢了不说还容易把更多噪声细节纳入模型中带来大量虚假极点。我的一般经验是三个原则i_blocks2采样时长要覆盖至少20个最低关注周期的数据i_blocks取数据总采样点数的5%到10%之间通常10到30之间足够在数据充裕的情况下尽量多试几组取值看结果是否稳定对上面的仿真数据总采样点数12000i_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对比计算效率、或者引入自动聚类算法让稳定图全自动判读都是价值很高的延伸方向。