做结构模态测试的人手里如果已经有几组加速度响应数据又不想被频域方法的各种窗函数和平均次数搞得心烦那么SSI-COV协方差驱动随机子空间识别是一个非常值得掌握的工具。它直接用环境激励下的响应数据来识别模态频率、阻尼比和振型不需要人工激励也不需要知道激励的具体大小。这篇文章就从工程应用的角度讲讲这套方法的原理、Matlab实现细节、以及我实际调参数踩过的坑。1. 背景为什么选SSI-COV而不是频域方法先聊一个很多人问过的问题都有了频响函数和峰值拾取法为什么还要用SSI-COV原因在于实际结构测试里激励往往不是可控的例如桥梁上的车辆荷载、建筑物上的风荷载、机械运行时的环境振动这些都是随机激励输入无法准确测量。频域方法需要把激励和响应做互谱或频响函数但环境激励下很难拿到高质量的激励信号强行用峰值法只能看到共振峰阻尼比算出来误差也大振型还可能被密集模态污染。SSI-COV走的是另一条路直接用响应的协方差序列来构造系统状态空间模型再把模型转化成模态参数。它的名字里“协方差驱动”指的就是这一步。它的好处非常实际不需要已知输入环境激励即可利用随机激励下输出响应的统计特性能同时识别频率、阻尼比、振型并且不像FDD那样靠峰值猜测模态阶次抗噪声干扰能力较强尤其对低频、弱模态比频域法稳定有成熟的状态空间理论基础模态参数的不确定性可以做后验评估当然它也有代价。最直观的就是计算量比频域方法高而且需要人为确定系统阶次。阶次选不好就很容易出现虚假模态。但这些问题通过稳定性图可以部分缓解我在后文会专门讲。我在参与的一个人行天桥振动测试项目里实测数据只有12个通道的加速度时程采样率256Hz总时长20分钟。用PolyMAX和SSI-COV都跑了一遍SSI-COV识别出的前四阶阻尼比虽然只有0.6%到2.3%但重复测试间的方差明显比频域方法小尤其第二阶和第三阶两个频率相距只有3Hz左右的模态SSI-COV能把振型分开峰值法则几乎糊成一团。那个项目之后我就把SSI-COV放进了常规工具箱。2. 核心原理状态空间模型和协方差序列的来龙去脉SSI-COV的理论基础是线性时不变系统的状态空间描述。一个多自由度结构在离散时间下的状态方程和观测方程可以写成x(k1) A * x(k) w(k) y(k) C * x(k) v(k)这里的x是状态向量包含位移和速度或者对应的离散状态变量y是实测输出A是系统矩阵C是观测矩阵w和v分别是过程噪声和测量噪声。我们做模态识别的目标就是从实测输出y(k)中估计系统矩阵A的特征参数再从中解出振动系统的固有频率、阻尼比和振型。这里有个关键点随机子空间的假设w和v都是零均值白噪声并且与系统状态不相关。这是一个很强的假设但工程上大多数环境振动激励逼近这个条件。如果数据里有强谐波干扰或者明显的非平稳漂移请看第5节预处理的方法。协方差驱动的思路是这样的先定义输出协方差序列R(i) E[y(ki) * y(k)^T]然后用这些协方差块来构造一个分块Toeplitz矩阵。为什么可以这么做因为线性系统本身有马尔可夫性未来输出和过去输入的统计相关性中包含了系统动态的全部信息。具体推导过程不展开最终结论是Toeplitz矩阵可以分解为可观矩阵和可控矩阵的乘积也就是T O * G其中O是观测矩阵和系统矩阵构成的可观矩阵G是可控矩阵。对T做奇异值分解SVD就能得到系统矩阵A的估计。得到A之后对它做特征值分解复特征值对应系统的极点利用极点和采样时间可以换算频率和阻尼比而观测矩阵C乘以特征向量就得到振型。我知道很多第一次看到这部分的人会被这些矩阵绕晕但用一句话总结就是协方差序列把激励的影响“平均”掉了剩下的就是系统本身的状态演化规律。SVD把噪声子空间和信号子空间分开取信号子空间就可以提取模态参数。这就是SSI-COV的全部骨架其他都是对这个骨架的完善和数值优化。3. 数据类型和预处理直接决定识别质量的三个步骤很多人一开始做SSI-COV识别出来一塌糊涂第一反应是算法不行其实八成是数据预处理没做好。数据预处理对SSI-COV的影响比对频域方法的影响还要大。因为协方差序列对趋势项、直流偏移、高频噪声都非常敏感任何一个不当处理都会在后续的SVD中变成虚假模态。我一般按下面三步走。第一步是去趋势和去均值。加速度传感器如果有轻微的温漂或零漂数据里就会有一个时变趋势这个趋势在计算协方差时相当于一个低频强分量会在低阶模态附近造成假峰。建议先用多项式拟合并去除趋势项同时减去均值。但多项式阶数不用太高通常1到2阶就够了太高反而会吃掉真实的低频模态。第二步是低通滤波。原始数据里高频噪声越少需要的系统阶次就越低稳定图上的假模态也越少。不过滤波会改变相位对频域方法影响很大但SSI-COV是基于统计特性的对相位不是特别敏感所以可以在滤波后直接使用。关键是滤波器要选择零相位版本例如Matlab的filtfilt函数可以避免相位偏移。截止频率一般取你关心的最高分析频率的1.5倍左右例如关心到50Hz就滤到75Hz再高只会增加噪声。第三步是重采样和降采样。这一步常被忽略。SSI-COV的计算量和采样点数直接相关长数据高采样率会让汉克尔矩阵和Toeplitz矩阵的维度爆炸没必要。重采样可以在低通滤波后进行采样率降到分析最高频率的2.5到3倍即可。例如原采样率1024Hz关心最高频率30Hz完全可以直接降采样到100Hz计算量减少90%识别结果基本不变。我经常用resample函数它会自动做抗混叠滤波比我手动滤波后再抽取靠谱。预处理之后最好把数据做一下可视化检查看所有通道的时程量级是否在一个数量级。如果某个通道的均方根值比其他通道小几十倍先查传感器不是每个通道都强制参与识别。在SSI-COV里通道量级差异过大会让SVD的奇异值偏向大通道小通道的振型信息容易被淹没。必要时对每个通道做归一化但振型的绝对值就会丢失只能从相对值角度使用所以不是万不得已我不做归一化而是通过传感器标定和放大倍数设置让所有通道基本一致。4. Matlab代码实现从协方差矩阵到频率阻尼振型下面给出我在工程中使用的SSI-COV主程序框架代码通俗易懂核心步骤都加了注释。4.1 主流程代码框架function [fn, zeta, phi] ssi_cov(y, fs, ncols, nrows, order) % y : 输出响应矩阵每一列为一个测点通道 % fs : 采样频率 % ncols : Toeplitz矩阵的列分块数 % nrows : Toeplitz矩阵的行分块数 % order : 系统阶次偶数 [N, nch] size(y); % 计算协方差序列 R(i)i 0, 1, ..., nrowsncols-2 maxlag nrows ncols - 2; y_mean mean(y, 1); y y - y_mean; R zeros(nch, nch, maxlag1); for k 0:maxlag % 对应样本协方差滞后k a y(1:N-k, :); b y(1k:N, :); R(:,:,k1) (a * b) / (N-k); end % 构造分块Toeplitz矩阵 T zeros(nch*nrows, nch*ncols); for i 1:nrows for j 1:ncols lag ncols - j i - 1; % 需要仔细推导的索引 T((i-1)*nch1:i*nch, (j-1)*nch1:j*nch) R(:,:,lag1); end end % SVD分解 [U, S, V] svd(T, econ); % 根据系统阶次order截断 order2 order / 2; U1 U(:, 1:order); S1 S(1:order, 1:order); V1 V(:, 1:order); % 计算状态矩阵A的估计 O U1 * sqrt(S1); % 可观测矩阵估计 O_bar O(1:end-nch, :); % 删掉最后nch行 O_up O(nch1:end, :); % 删掉最前nch行 A_est O_bar \ O_up; % 观测矩阵C估计取可观测矩阵的前nch行 C_est O(1:nch, :); % 特征值分解 [Psi, Lambda] eig(A_est); lambda diag(Lambda); % 离散极点转连续极点 mu log(lambda) * fs; % 频率和阻尼比 fn abs(mu) / (2*pi); zeta -real(mu) ./ abs(mu) * 100; % 百分比% % 振型 phi C_est * Psi; % 每个特征向量的列对应一个模态 end这段代码可以跑通但我必须提醒几点。4.2 关于索引和矩阵维度的几个大坑第一Toeplitz矩阵的索引是最容易写错的。我上面的代码用的是lag ncols - j i - 1这个公式取决于R(1)对应滞后0。一旦把滞后次序搞反可能识别出的频率是正确的但阻尼比符号是反的振型也会乱。我自己的建议是不要凭记忆写这段用一个滞后矩阵从左到右、从上到下打出来检查一遍。第二系统阶次order必须为偶数因为每一阶物理模态对应一对共轭复极点。如果是奇数特征值分解后会出现一个实数极点对应的“模态”不是物理意义的共振后续很容易被当成虚假模态误删。第三A_est的求解用O_bar \ O_upMatlab默认的最小二乘。这一步如果O的条件数很差结果就会不稳。此时可以尝试在SVD截断前进一步做奇异值截断。S值如果出现跳崖式下降说明截断阶次应该选在跳崖点前的平坦区域。下面我专门讲怎么定阶。5. 系统阶次怎么定奇异值曲线和稳定图实战SSI-COV最让人纠结的就是order取值。order定太小模态会漏掉阻尼比估计偏差也大order定太大稳定图上全是计算产生的数值极点筛选起来头痛。我的做法是两层筛选结合。第一层看奇异值曲线。SVD得到的奇异值S对角元从大到小排列信号子空间的奇异值远大于噪声子空间的奇异值所以在对数坐标下奇异值曲线从陡降变成平缓的转折点就是有效阶次的参考位置。转折点之前的数量乘以2大致对应系统阶次上限。实际使用时我在那个点附近再扩大一倍留出余量比如转折点在25就可以尝试order在30到50之间变化。第二层用稳定图。Stabilization Chart我在实际项目里必做步骤很简单设定一个order序列比如从10到80步长为2对每个order跑一次SSI-COV得到一组频率、阻尼比、振型按相似容差判断哪些模态在多次识别中稳定稳定性的标准一般这样设置% 频率容差 1% f_tol 0.01; % 阻尼比容差 10%阻尼本身识别难度大可放宽 d_tol 0.10; % MAC值容差 2% mac_tol 0.02;如果某个模态候选在相邻两次order识别中频率变化小于1%阻尼比变化小于10%相对值MAC大于98%就可以认定为稳定点。把所有稳定点画在“频率-order”平面上形成竖线竖线聚集处就是真实模态。稳定图看起来直观但筛得很科学。值得注意的是阻尼比容差不能太紧。阻尼比本身是识别量里噪声最大的一个尤其环境激励下2%的阻尼比和2.2%的阻尼比在物理上可能没区别但相对变化超过10%容易把稳定点漏掉。我通常把阻尼的稳定性判据放宽到20%频率还是1%。下面给一个画稳定图的常见代码段orders 10:2:60; all_fn cell(length(orders), 1); all_zeta cell(length(orders), 1); all_phi cell(length(orders), 1); for oi 1:length(orders) [fn_o, zeta_o, phi_o] ssi_cov(y, fs, 20, 20, orders(oi)); % 只保留正频率且阻尼比在0到10%之间的假想模态 valid fn_o 0 zeta_o 0 zeta_o 10; all_fn{oi} fn_o(valid); all_zeta{oi} zeta_o(valid); all_phi{oi} phi_o(:, valid); end % 循环配对并标记稳定点 % 具体实现可根据数据量优化这里只给出思虑框架稳定图的绘制代码不难但耗时在配对逻辑上。如果每个order识别出20个候选模态60个order就是1200次配对用for循环也能跑但数据量大时要向量化。我在代码里会先用频率排序做预筛再用MAC矩阵一次算完这样能省很多时间。6. 实测踩坑记录为什么识别的阻尼老是漂阻尼比对SSI-COV来说是个难点。我整理了自己项目里常见的问题列成一张表大家可以直接对照排查。现象可能原因排查和处理高频模态识别出很多假极点系统阶次过高、未充分滤波降低order上限提高截止频率的滤波质量低频段出现负阻尼趋势项去除不彻底或数据截断效应检查去趋势增加重采样后的数据长度各阶阻尼比几乎都一样存在强噪声通道或某个通道饱和检查时程剔除异常通道后再识振型向量相位断续乱跳传感器方向接反或符号定义不一致检查传感器方向统一坐标方向约定稳定图上竖线很多但MAC不高噪声过大或测点布置不敏感增加平均次数换句话说增加数据长度还有一个我一开始没注意的问题是数据长度对阻尼识别影响极大。SSI-COV对协方差的估计质量依赖于数据长度。理论上协方差序列R(i)需要足够多的样本才能收敛。如果数据长度只有几十秒低频模态的协方差还没被充分平均阻尼比的估计会偏向于随机游走。我现在的经验是最低确保每个分析频带内至少500个循环周期。比如关注最低频率1Hz至少采集500秒如果最低频率0.1Hz至少采集5000秒。这在实际桥梁测试中往往要熬夜挂机。另一个容易出错的是重采样导致的虚假模态。我遇到过识别出一个12.5Hz的模态怎么调参数都不消失后来发现是原始数据里面有50Hz市电工频干扰resample到100Hz后混叠到了12.5Hz附近。解决方案是先做50Hz陷波滤波或者先用零相位低通滤到远低于奈奎斯特频率后再降采样。这个例子充分说明预处理做的多细致后续就能多省心。另外如果要识别多参考点振型并希望振型归一化到某一点建议在输入SSI-COV之前就记录好测点坐标和通道对应关系。识别之后振型相位是复数域的信息土木结构小阻尼情况下相位接近0°或180°但如果识别出的复数振型相位在空间上不是连续分布往往不是结构问题而是传感器方向标错了。我在一次试验中吃过亏后来养成了识别完成后画振型动画的习惯比单纯看数字可靠得多。7. 自动筛选真实模态MAC和模态参与因子的使用稳定图帮我们缩小了范围但真实工程数据里依然会有一些“看起来稳定”但物理上说不通的极点比如数值模态和真实模态频率太接近或者谐波分量近似周期信号造成的伪极点。这时候需要另外两个工具MACModal Assurance Criterion和模态参与因子。MAC是判断两个振型向量相关程度的工具公式是MAC_ij abs(phi_i * phi_j)^2 / ((phi_i*phi_i) * (phi_j*phi_j));当两个振型来自不同结构模态时MAC值理论上接近0数值上一般小于0.2来自同一模态时接近1。在自动筛选时我会这样做计算所有候选模态之间的MAC矩阵把MAC超过0.9且频率也接近的模态合并为一类每一类里选稳定图上出现次数最多的那个频率作为最终模态这个步骤可以有效去除重复识别。但MAC也有局限如果传感器数量少或者测点分布不理想不同振型之间的MAC可能天然很高。比如只有两个测点且对称布置反对称模态和对称模态在测点上就可能分不开。此时只能靠频率差和阻尼比经验判断没法完全自动化。模态参与因子主要用来识别谐波分量。环境激励下的旋转机械或桥梁吊杆经常会产生近似正弦的窄带激励。这种激励在SSI-COV里会使某个极点的“参与因子”特别大而且这个极点对应频率的奇异值也很大。参与因子可以从可控矩阵的估计值里算出来但由于代码篇幅原因不展开。实际项目中我会在目标频率附近留出一段“警戒区”把任何频率落在电网工频、转速基频附近的候选模态都标记为待验证配合现场工况记录判断是否剔除。8. 实用工具箱推荐和代码调优思路如果不想从零开始写SSI-COVMatlab里有几个现成的工具箱可供参考我在项目中也混用过MACEC比利时鲁汶大学开发的结构模态分析工具箱里面有SSI-COV和稳定图实现学术圈用得很多可以参考源码学习细节。MATLAB System Identification Toolboxn4sid和ssest函数可以实现部分子空间辨识但面向一般线性系统输出的是状态空间模型需要自己做模态解算。一些GitHub上的开源实现搜索cov-ssi或stochastic-subspace质量参差不齐用之前一定要检查Toeplitz索引和奇异值截断逻辑。我的建议是第一次做这个研究的人先自己写一个简化版比如只识别单自由度系统的模态参数把每个矩阵打印出来逐步理解。这样比直接拿工具箱跑真实数据学习效率高得多。理解原理之后再考虑封装成自己的函数库因为实际项目的通道数、采样频率、阶次范围差异很大开源代码不一定能直接复用。调优时有几个思路值得尝试分频带处理。如果关心0.1Hz到100Hz很宽的频带一次做完往往效果不好。可以先低通滤到30Hz识别低频段再把高通滤掉低频部分识别高频段但要注意滤波造成的边界效应数据首尾要预留足够长度的丢弃区。多段数据平均。如果现场只记录了两次工况可以把两次数据拼接后进行整体识别吗不建议直接拼接因为两次数据的环境激励水平和结构状态可能不同。更好的做法是分别识别再将频率和振型结果做加权平均权重和信噪比相关。通道加权。如果测点较多某些通道信噪比明显偏低可以在构建协方差矩阵之前对数据做加权但加权会影响振型的绝对值因此我只在严重不平衡时使用而且识别完会把振型重新缩放到实际测点单位。工具箱不是重点重点是你能不能用里面的函数解释清楚每一个输入参数的含义。曾经有朋友问我为什么他用了某个现成代码后阻尼比每天跑结果都不一样我让他看了代码里对协方差滞后数的选择。那一版代码里ncols和nrows写死了根本不随数据长度变化低频部分协方差估计方差极大阻尼比自然天天变。这种问题只有自己掌握了原理才看得懂。9. 一次真实数据的完整识别流程复盘我用一个人行天桥案的例子把整套流程串起来讲方便大家照葫芦画瓢。那次测试一共布了8个加速度传感器分布在主梁第二跨的八等分点上采样率200Hz连续记录20分钟。我按下面的流程走了一遍。数据预处理fs 200; % 原始数据矩阵 y_raw每列一个通道 y_raw detrend(y_raw, 1); % 线性去趋势 y_raw y_raw - mean(y_raw, 1); % 设计零相位低通截止40Hz [b, a] butter(4, 40/(fs/2), low); y_fit filtfilt(b, a, y_raw); % 降采样到100Hz fs_new 100; y_new resample(y_fit, fs_new, fs);参数选择上nrows和ncols我取20。注意这个参数对应Toeplitz矩阵里用到的协方差滞后数不是越高越好。理论上滞后阶数越大包含的模态信息越多但协方差估计误差也随滞后增大而增大。常看到的建议是ncols取最大滞后对应10到30个采样周期覆盖最关心模态的3个周期以上。以最低模态频率1Hz为例100Hz采样下对应100个采样点为1周期ncols20对应滞后到2s足够覆盖。但如果是0.2Hz的超低频结构ncols20只覆盖1s远远不够必须加大。这个关系很好理解。之后对order从20到80步长2跑了稳定图发现三列很明确的稳定线分别对应约1.9Hz、4.6Hz、8.2Hz。第四阶8.2Hz在较大order时才稳定阻尼比约1.1%。最终以order42的频率结果作为最终输出因为该阶次下所有目标模态都稳定且没有多余假极点。振型结果输出前我把8个测点的相对振幅按测点坐标做了插值画成侧视图确认第二阶振型在中跨出现反弯点和有限元预测一致。这一步不可省因为只有通过振型空间形态的物理合理性检查才能确认识别结果不是数值伪模态。阻尼比的结果在三次重复时间段里分别是一阶0.8%、0.9%、0.8%二阶1.5%、1.6%、1.4%三阶1.0%、1.1%、1.0%试想如果不做多段重复一次测试的阻尼比随机误差可能超过20%。这类统计稳定性信息在结构健康监测里的阈值设置和模型修正中非常重要。10. 后续还能怎么扩展SSI-COV不是一个只能离线处理的算法现在已经有不少在线识别方案用滑动窗口实时更新协方差矩阵和模态参数也有将它与贝叶斯方法结合输出模态参数的概率分布还有人在做自动化OoMAOperational Modal Analysis把稳定图自动生成、自动筛选、自动报告封装成一套流程。对于Matlab使用者来说从本文框架出发扩展并不难把Toeplitz矩阵构造改成稀疏存储可以支持更大的数据量和更多通道把SVD替换为随机化SVD计算速度能提升一个量级适合无人机巡检采集的大量数据与有限元模型结合通过MAC和频率残差做模型修正这是结构健康监测领域很热的方向我个人建议如果研究方向是结构健康监测SSI-COV这一套方法论不仅要会调用更要亲手推导一遍状态空间模型和协方差Toeplitz分解。这个过程能把很多零散知识串起来也不至于在数据异常时完全摸不着头脑。最后分享一个小体会识别结果出来后先别急着相信计算机给的那几个数。我在实际项目中已经养成习惯把识别出的频率和理论估算先对一眼把振型形状画出来和结构几何形状对比一下。凡是和物理直觉冲突很大的结果九成是输入数据或参数设置出了问题而不是结构真的出了什么奇异的模态。用这套思路搭配SSI-COV做起模态分析来会稳很多。
