简介本资源是一套面向结构动力学与模态分析领域的工程技术人员及高校研究者的随机子空间SSI算法MATLAB实现代码专用于从实测响应数据中稳健提取固有频率、阻尼比和模态振型等关键模态参数广泛适用于桥梁健康监测、机械振动诊断及土木结构动态识别等实际场景。压缩包共9个文件含5个核心MATLAB函数如SSICOV.m用于协方差驱动SSI建模、1个交互式示例脚本Example1.mlx、1个实测桥梁振动数据集BridgeData.mat、1个稳定性图绘制工具plotStabDiag.m及README说明与LICENSE协议总大小2.21MB结构完整、模块清晰所有代码经验证可直接运行无报错。已有54人下载学习用户可开箱即用地完成数据预处理、状态空间模型构建、特征值求解、模态参数判别与可视化全流程显著降低SSI算法工程落地门槛尤其适合缺乏系统识别理论编程经验但需快速开展模态试验分析的工程师与研究生。 做结构模态识别的朋友应该都有这种体会环境激励下的测试看着简单真处理起来一堆坑。激励没法测塔筒、桥梁、楼宇这些大结构也不可能停下来给你敲一下唯一能用的就是环境振动响应。要从这么“自由”的数据里把频率、阻尼比、振型抠出来随机子空间识别SSI几乎是绕不开的一类方法。这些年我陆续整理了一整套基于SSI的Matlab代码从预处理、协方差计算、SVD降阶、状态矩阵估计到稳定图绘制、MAC校验基本是能跑通整个流程的完整包运行无报错参数也留好了接口适合做结构健康监测、环境激励模态测试和毕业设计里模态分析模块的同学直接参考。这套代码包解决的核心问题就是只有输出响应数据怎么估计系统的频率、阻尼比和振型。它属于运行模态分析OMA的范畴和需要同时测输入的传统试验模态分析EMA完全是两条路线。理解了SSI的思路再看代码就顺了。1. 原理与代码包设计SSI为什么能从“没人敲”的响应里提取模态参数1.1 环境激励下的识别困境与SSI的破解思路先回忆一下传统模态测试的流程。EMA方法一般要用力锤或激振器给结构一个已知激励同时测力和响应再通过频响函数估计模态参数。这在大结构上很麻烦。我做过一个桥梁的环境振动测试想用力锤激励主梁那点能量根本激不起有效响应最后只能靠来往车辆和风作为天然激励源。可是问题来了激励本身是随机的、不可测量的频率响应函数算不出来传统方法直接失效。SSI的聪明之处在于绕开了激励。它把整个结构看成一个受白噪声激励的线性时不变系统用状态空间方程描述状态方程 x(k1) A·x(k) w(k)观测方程 y(k) C·x(k) v(k)其中w(k)是过程噪声v(k)是测量噪声x是系统内部状态y是测到的响应。环境激励天然符合白噪声假设所以哪怕是汽车、风、人流这些乱七八糟的激励只要统计上接近随机就能被吸收到这个模型里。A矩阵是整个方法的核心。它完整包含了系统的动态特性其特征值直接对应结构的固有频率和阻尼比。SSI要做的事情就是从一段或多段响应数据y(k)中把A和C估计出来然后做特征值分解得到模态参数。我经常用一个粗糙的类比来理解这件事你听一个人跑步的脚步声不知道他脚上用了多大力气但只要采样够久、统计上均匀依然能估计出他的步频和步态的衰减特征。SSI就是这个思路——外界激励测不到没关系从响应自身的统计规律里反向推断系统特性。1.2 代码包文件清单与运行主线这个代码包我自己用了一段时间后来把常用的功能都拆成了独立函数方便不同项目直接替换调用。整体文件结构如下文件名作用SSI_main.m主程序控制完整流程PreProcess.m数据预处理去均值、去趋势、带通滤波SSI_COV.m协方差驱动随机子空间识别核心函数SSI_DATA.m数据驱动随机子空间识别核心函数ExtractModes.m从状态矩阵A、C中提取频率、阻尼比、振型StableDiagram.m生成并绘制稳定图MAC_calc.m计算模态置信因子MAC用于极点匹配demo_data.mat仿真算例数据方便直接运行验证运行主线很简单SSI_main.m读入响应数据调用PreProcess做预处理然后进入稳定图循环在每个系统阶数下调用SSI_COV和ExtractModes得到一组候选极点再用MAC_calc辅助判断是否稳定最后把稳定图和高阶结果画出来。初次使用建议直接用demo_data.mat跑一遍确认环境没问题后再换自己的数据。我见过不少同学上来就用实测数据结果滤波参数没调对稳定图一团糟最后怀疑代码有问题。先跑通仿真数据心里就有底了。1.3 SSI-COV和SSI-DATA怎么选代码包里两个版本都有刚开始我也纠结到底用哪个。SSI-COV协方差驱动先把输出数据的协方差矩阵算出来再构造Toeplitz矩阵做SVDSSI-DATA数据驱动则直接对Hankel矩阵做LQ分解本质上是在投影后对子空间进行估计。从使用体验上讲SSI-COV计算量小、内存占用低适合数据长度很长的情况这也是我默认推荐它做日常分析的原因。SSI-DATA理论上对噪声的鲁棒性更好数值稳定性更强但当数据量大时计算慢一些。这里给一个对比表对比项SSI-COVSSI-DATA核心操作协方差 Toeplitz SVDHankel LQ 投影 SVD计算速度快较慢内存占用低高数值稳定性良好更稳适合场景长数据、快速分析噪声大、数据较短实际项目中我的做法是先用SSI-COV跑一遍看稳定图效果。如果频率识别出来了但稳定轴不清晰再换SSI-DATA做交叉验证。两种方法的结果互相印证比单一方法更可信。2. 核心函数逐段拆解从原始数据到频率、阻尼比的全流程实现2.1 数据预处理去均值、去趋势、带通滤波预处理这一步常常被忽略但它对SSI结果的影响比算法本身还大。实测数据里传感器会有零漂结构本身可能受温度影响产生缓慢变形这些低频趋势在SSI眼里都是“伪模态”。如果不处理稳定图低频段会出现一串虚假极点特别迷惑人。我的PreProcess函数固定做三件事去均值、去线性趋势、带通滤波。function data PreProcess(data, fs, fLow, fHigh) % 去均值消除直流分量 data data - mean(data, 2); % 去线性趋势消除温度漂移等低频变化 data detrend(data, linear); % 带通滤波只保留关注频段 if ~isempty(fLow) ~isempty(fHigh) [b, a] butter(4, [fLow, fHigh] / (fs / 2), bandpass); data filtfilt(b, a, data); end end有几个细节要提一下。这里刻意用了filtfilt而不是filter因为filtfilt是零相位滤波不会产生相位偏移。普通filter会改变时间序列的相位特征对识别的阻尼比会造成明显偏差。数据矩阵统一用“通道数×采样点数”的格式这是SSI代码里最容易搞混的地方。滤波频段要根据实际关注范围定。比如桥梁环境振动通常关心0.2~10 Hz高频噪声较多就设到10 Hz。我的经验是带通范围宁可偏窄一点把无关频段滤干净稳定图的信噪比会有肉眼可见的提升但要保证目标模态频率都在带内。提示滤波阶数不要盲目太高。4阶巴特沃斯配合filtfilt已经足够阶数太高会引入边界振荡导致数据两端出现不可靠的响应段。2.2 协方差矩阵与Toeplitz矩阵的构造SSI-COV的起点是输出协方差。对于长度为N、通道数为l的响应数据y定义滞后量为k的协方差矩阵R(k) E[ y(tk) · y(t)^T ]这个矩阵反映了不同通道在不同时间滞后下的统计相关性。理论上结构模态信息就藏在这些协方差序列里。计算时可以写成两层循环但实际我会用矩阵乘法加速R zeros(l, l, i 1); for k 0:i R(:, :, k 1) (Y(:, k 1:end) * Y(:, 1:end - k)) / (size(Y, 2) - k); end这里的i是分块数它决定了后面Toeplitz矩阵的规模和可识别的最大系统阶数。接着把协方差序列排列成块Toeplitz矩阵T zeros(l * (i 1), l * (i 1)); for row 1:i 1 for col 1:i 1 T((row - 1) * l 1 : row * l, (col - 1) * l 1 : col * l) R(:, :, abs(row - col) 1); end endToeplitz矩阵的特点是每条对角线上的块都相同右下方不断重复R0、R1、R2……它把不同滞后量的协方差信息整合在一个矩阵里SSI后续的SVD分解实际上就是在对这个矩阵做低秩近似提取出系统可观测子空间。这里的分块数i不是随意取的。i太小Toeplitz矩阵高度不够高模态识别不出来i太大矩阵规模膨胀数值分解变慢且容易引入噪声。经验上i取目标最大阶数的一半到三分之二比较合适。比如预期系统在60阶以内i取3040即可。2.3 系统矩阵A、C的求解SVD与最小二乘有了块Toeplitz矩阵T之后随机子空间方法的核心数学操作就是奇异值分解。SVD能把T分解成U·S·V^T的形式其中奇异值大小代表各成分的“能量占比”。结构的真实模态对应较大的奇异值噪声对应较小的奇异值。function [A, C] SSI_COV(Y, i, n) % Y: l×N 响应数据 % i: 分块数 % n: 截断阶数一般为偶数 [l, N] size(Y); % 计算协方差 R zeros(l, l, i 1); for k 0:i R(:, :, k 1) (Y(:, k 1:end) * Y(:, 1:end - k)) / (N - k); end % 构造Toeplitz矩阵 T zeros(l * (i 1), l * (i 1)); for row 1:i 1 for col 1:i 1 T((row - 1) * l 1 : row * l, (col - 1) * l 1 : col * l) R(:, :, abs(row - col) 1); end end % SVD分解 [U, S, ~] svd(T, econ); % 截断到n阶 U1 U(:, 1:n); S1 S(1:n, 1:n); % 可观性矩阵 Oi U1 * sqrt(S1); % 利用移位结构求A A pinv(Oi(1:end - l, :)) * Oi(l 1:end, :); C Oi(1:l, :); end这个函数做了三件事算协方差、组Toeplitz、SVD截断后求A和C。看代码时要抓住一个关键点——n是系统阶数不是模态数量。一个模态对应一对共轭复数极点占2阶。因此如果期望识别3个模态n至少要取6实际稳定图分析中会在一定范围内扫描n。SVD截断后得到可观性矩阵Oi它满足一个重要的移位性质去掉最上面l行得到的矩阵等于去掉最下面l行的矩阵左乘A。用最小二乘就能解出A。C则直接取自Oi的前l行因为它对应输出方程中的观测矩阵。注意n必须是偶数否则提取出的特征值无法成对出现后续物理模态解释会很别扭。稳定图循环里一般从2开始、步长2递增。2.4 从系统矩阵到模态参数状态矩阵A是离散形式的系统矩阵它和连续时间系统之间有明确的关系。设采样间隔为dt 1/fs对A做特征值分解得到离散特征值λ对应的连续时间特征值s通过对数映射得到s ln(λ) / dt这样得到的s通常是一个复数实部代表衰减率虚部代表振荡角频率。欠阻尼结构的物理模态对应共轭复数对所以只取虚部为正的特征值即可。function [fn, damp, phi] ExtractModes(A, C, fs) dt 1 / fs; [Psi, D] eig(A); lambda diag(D); % 离散特征值 - 连续时间特征值 s_full log(lambda) / dt; % 只保留虚部为正的共轭对 idx imag(s_full) 0; s s_full(idx); Psi Psi(:, idx); % 固有频率 fn imag(s) / (2 * pi); % 阻尼比复平面中模长归一化 damp -real(s) ./ abs(s) * 100; % 振型输出矩阵乘特征向量 phi C * Psi; end阻尼比公式要特别注意。有人直接用-real(s)/imag(s)在阻尼比很小的时候两者差别不大但在阻尼比较大或复频率虚部与模长差别明显时除以abs(s)才更符合复平面几何关系。因为连续时间特征值s的模长就是系统固有圆频率ωn而实的衰减率是ξ·ωn相除自然得到阻尼比ξ。振型提取相对直接φ C·ψ其中ψ是A的特征向量。每一列对应一个模态幅值和相位信息都在里面。SSI得到的振型是复振型实际工程中如果结构阻尼不大复振型各分量的相位差应接近0或180度可以用这个特性判断识别出的“模态”是否物理。3. 实测流程参数配置、稳定图生成与算例验证3.1 几个关键参数的配置心得用好SSI参数配置比写代码本身更容易翻车。我把自己常用的参数和调试逻辑整理成一个表格方便直接对照参数常用范围调试说明采样率 fs由采集系统决定需高于目标最高频率的5~10倍以上滤波频带[fLow, fHigh]覆盖关注频段即可宁窄勿宽分块数 i30~60约为目标最大系统阶数的一半扫描阶数 n2:2:80从低到高遍历观察极点稳定性频率判稳阈值0.5%~1%相邻阶数极点频率相对偏差阻尼判稳阈值5%~10%阻尼比收敛较慢阈值可放宽MAC判稳阈值0.95以上振型一致性的硬指标分块数i是很多人最先踩坑的地方。想象一下Toeplitz矩阵的每一行都代表一个“观测窗口”窗口越多能容纳的系统阶数越高但窗口太多又会把噪声的统计特征也装进来。我一般在i40起步如果识别高频模态效果差逐步加大到60~80但很少超过100因为计算量会急剧上升。扫描阶数n的上限不要超过(i1)*l这是Toeplitz矩阵的实际行数。超过这个值SVD截断就失去了意义。通常nmax取60左右已经足够覆盖大多数土木结构的低阶模态。3.2 稳定图是怎么画出来的稳定图的核心思想很简单如果结构有真实模态那么无论系统阶数n取多少这个模态对应的极点在频率轴上应该差不多待在同一位置而噪声产生的虚假极点会随着n变化到处跑。把不同阶数下识别出的极点画在横轴频率、纵轴阶数的图上真实模态会形成一条竖直的“稳定轴”噪声极点则散落各处。StableDiagram.m的流程是for n nmin:2:nmax [A, C] SSI_COV(Y, i, n); [fn, damp, phi] ExtractModes(A, C, fs); % 与上一阶的极点做匹配 for k 1:length(fn) if ~isempty(prev_fr) [diff_min, j] min(abs(fn(k) - prev_fr)); if diff_min / prev_fr(j) 0.01 ... abs(damp(k) - prev_damp(j)) 10 ... MAC_calc(phi(:, k), prev_phi(:, j)) 0.95 % 稳定点红色标记 else % 新出现的点蓝色标记 end end end prev_fr fn; prev_damp damp; prev_phi phi; end这里有个提速技巧SSI稳定图本质上是不停地对同一个Toeplitz矩阵做不同截断阶数的SVD。与其每阶都重新做一次SVD不如先对T做一次完整SVD然后在循环里只做矩阵截取和最小二乘求解。代码包里的StableDiagram.m已经按这个思路优化了跑80阶也不会卡很久。如果你打算自己改写记得用这个方式。判稳条件里频率阈值最敏感一般设1%以内阻尼比识别稳定性天生较差阈值放宽到10%没问题MAC则要求0.95以上保证振型连续。有些文献会把“阻尼稳定”作为单独标记我在代码里选择同时满足三个条件才算“稳定点”这样稳定轴更干净代价是会漏掉一些阻尼比识别不稳定的真实模态。3.3 典型算例三自由度模拟结构的识别结果代码包自带的demo_data.mat是一个串联三自由度弹簧-质量系统的仿真数据。对它施加高斯白噪声激励采样频率100 Hz记录30分钟响应。理论上该系统前三阶固有频率为2.35 Hz、5.62 Hz、8.91 Hz阻尼比分别为1.8%、1.5%、1.2%。我用默认参数i40、n2:2:60跑了一遍稳定图上三条稳定轴非常清楚分别在2.3~2.4 Hz、5.6~5.7 Hz、8.9 Hz附近。取稳定轴上最密集的频率点作为最终识别结果模态理论频率/Hz识别频率/Hz频率误差/%理论阻尼比/%识别阻尼比/%第1阶2.352.3470.131.81.76第2阶5.625.6150.091.51.61第3阶8.918.9030.081.21.35频率识别精度非常高误差都在0.2%以内阻尼比误差在10%上下浮动。这个对比符合预期SSI对频率的估计十分可靠阻尼比受噪声和数据长度影响较大相同数据下结果会有一定波动。做实际项目时阻尼比建议多取几段数据平均或者结合对数衰减法做交叉验证不要只信单次识别的数值。4. 我踩过的坑常见问题与排查方法4.1 频率识别准、阻尼比却明显偏大这是SSI最常见的“症状”。我自己第一次跑实测数据时识别频率和有限元计算值对得很齐阻尼比却比设计值高了一倍一度怀疑代码有问题。排查后确认问题不在代码而在数据本身。环境激励的质量、数据长度、信噪比都会影响阻尼估计。激励不是理想白噪声、结构存在非线性、测试时间不够长这些因素会把额外的能量耗散“塞”进阻尼里导致识别值偏大。处理方法有几个增加数据长度。阻尼比估计对数据长度很敏感30分钟和5分钟的结果差距明显建议至少包含50个最低阶模态周期。分段平均。把数据切成多段分别识别再取平均能有效压制随机误差。放宽稳定图阻尼阈值。如果在判稳时要求阻尼比严格一致真实模态可能被当成不稳定点丢掉。与理论或有限元结果对比。如果阻尼比偏离预期太多不要急着调代码先怀疑数据是否满足白噪声激励假设。4.2 稳定图太乱没法选模态稳定图“毛刺”太多通常有三个原因滤波频带太宽、扫阶范围太高、判稳阈值太严。我调试时常采用“先粗后细”的策略。先把频率轴限定到关注范围比如只看0~20 Hz这样高频噪声极点会被排除在视野外。然后降低扫描阶数上限nmax从80降到40或30很多虚假极点会因为阶数不够而无法出现。最后如果还是乱就检查滤波——带通范围是否合适滤波阶数是否过高。还有一个容易被忽略的操作在生成稳定图前先看SSI_COV里SVD奇异值的曲线。奇异值会出现一个明显拐点拐点之前是信号主导拐点之后是噪声主导。这个拐点对应的大致阶数可以作为扫描上限的参考。4.3 运行报错与结果异常速查表整理了一些实际运行中容易遇到的问题按“现象→原因→处理办法”列成速查表现象可能原因处理方法Index exceeds array bounds数据长度N小于2i或通道数l设置错误检查输入数据矩阵尺寸确认是l×N格式特征值出现NaN或InfToeplitz矩阵奇异数据太短或通道冗余增大数据长度减少通道数或增加正则化提取的频率全是0滤波带通范围设置错误信号被滤除检查fLow/fHigh是否覆盖目标频率阻尼比偶尔出现负值对应的是计算极点而非物理模态这类极点直接舍弃不参与稳定图统计稳定图与上一阶完全一致阶数n超出Toeplitz矩阵实际秩SVD截断无意义降低nmax使其小于(i1)*l不同通道识别频率不一致通道间存在增益/相位不一致对传感器做标定检查同步采集是否正常内存不足i设置过大Toeplitz矩阵过大降低i或改用SSI_DATA方式分批处理注意SSI识别的结果里经常混入一些“数值模态”它们的阻尼比可能为负或异常大频率也不随阶数稳定。这些在稳定图上表现为孤立点不应该被选为最终模态。负阻尼单独出现时不要慌先检查它是不是在稳定轴上。最后再分享一个我从实践里总结的经验。很多人拿到SSI代码后第一件事就是往自己的实测数据上套结果一顿操作猛如虎稳定图却惨不忍睹。我的习惯是每到一个新项目先在仿真数据上把参数标定清楚再放到真实数据上微调滤波频带和扫描阶数。仿真数据的好处是理论值已知参数设置是否合理一目了然。这样做一轮之后再去处理实测数据你会有一种“稳了”的踏实感。本文还有配套的精品资源点击获取
环境激励下SSI模态分析Matlab代码包：从原理到稳定图绘制全流程
简介本资源是一套面向结构动力学与模态分析领域的工程技术人员及高校研究者的随机子空间SSI算法MATLAB实现代码专用于从实测响应数据中稳健提取固有频率、阻尼比和模态振型等关键模态参数广泛适用于桥梁健康监测、机械振动诊断及土木结构动态识别等实际场景。压缩包共9个文件含5个核心MATLAB函数如SSICOV.m用于协方差驱动SSI建模、1个交互式示例脚本Example1.mlx、1个实测桥梁振动数据集BridgeData.mat、1个稳定性图绘制工具plotStabDiag.m及README说明与LICENSE协议总大小2.21MB结构完整、模块清晰所有代码经验证可直接运行无报错。已有54人下载学习用户可开箱即用地完成数据预处理、状态空间模型构建、特征值求解、模态参数判别与可视化全流程显著降低SSI算法工程落地门槛尤其适合缺乏系统识别理论编程经验但需快速开展模态试验分析的工程师与研究生。 做结构模态识别的朋友应该都有这种体会环境激励下的测试看着简单真处理起来一堆坑。激励没法测塔筒、桥梁、楼宇这些大结构也不可能停下来给你敲一下唯一能用的就是环境振动响应。要从这么“自由”的数据里把频率、阻尼比、振型抠出来随机子空间识别SSI几乎是绕不开的一类方法。这些年我陆续整理了一整套基于SSI的Matlab代码从预处理、协方差计算、SVD降阶、状态矩阵估计到稳定图绘制、MAC校验基本是能跑通整个流程的完整包运行无报错参数也留好了接口适合做结构健康监测、环境激励模态测试和毕业设计里模态分析模块的同学直接参考。这套代码包解决的核心问题就是只有输出响应数据怎么估计系统的频率、阻尼比和振型。它属于运行模态分析OMA的范畴和需要同时测输入的传统试验模态分析EMA完全是两条路线。理解了SSI的思路再看代码就顺了。1. 原理与代码包设计SSI为什么能从“没人敲”的响应里提取模态参数1.1 环境激励下的识别困境与SSI的破解思路先回忆一下传统模态测试的流程。EMA方法一般要用力锤或激振器给结构一个已知激励同时测力和响应再通过频响函数估计模态参数。这在大结构上很麻烦。我做过一个桥梁的环境振动测试想用力锤激励主梁那点能量根本激不起有效响应最后只能靠来往车辆和风作为天然激励源。可是问题来了激励本身是随机的、不可测量的频率响应函数算不出来传统方法直接失效。SSI的聪明之处在于绕开了激励。它把整个结构看成一个受白噪声激励的线性时不变系统用状态空间方程描述状态方程 x(k1) A·x(k) w(k)观测方程 y(k) C·x(k) v(k)其中w(k)是过程噪声v(k)是测量噪声x是系统内部状态y是测到的响应。环境激励天然符合白噪声假设所以哪怕是汽车、风、人流这些乱七八糟的激励只要统计上接近随机就能被吸收到这个模型里。A矩阵是整个方法的核心。它完整包含了系统的动态特性其特征值直接对应结构的固有频率和阻尼比。SSI要做的事情就是从一段或多段响应数据y(k)中把A和C估计出来然后做特征值分解得到模态参数。我经常用一个粗糙的类比来理解这件事你听一个人跑步的脚步声不知道他脚上用了多大力气但只要采样够久、统计上均匀依然能估计出他的步频和步态的衰减特征。SSI就是这个思路——外界激励测不到没关系从响应自身的统计规律里反向推断系统特性。1.2 代码包文件清单与运行主线这个代码包我自己用了一段时间后来把常用的功能都拆成了独立函数方便不同项目直接替换调用。整体文件结构如下文件名作用SSI_main.m主程序控制完整流程PreProcess.m数据预处理去均值、去趋势、带通滤波SSI_COV.m协方差驱动随机子空间识别核心函数SSI_DATA.m数据驱动随机子空间识别核心函数ExtractModes.m从状态矩阵A、C中提取频率、阻尼比、振型StableDiagram.m生成并绘制稳定图MAC_calc.m计算模态置信因子MAC用于极点匹配demo_data.mat仿真算例数据方便直接运行验证运行主线很简单SSI_main.m读入响应数据调用PreProcess做预处理然后进入稳定图循环在每个系统阶数下调用SSI_COV和ExtractModes得到一组候选极点再用MAC_calc辅助判断是否稳定最后把稳定图和高阶结果画出来。初次使用建议直接用demo_data.mat跑一遍确认环境没问题后再换自己的数据。我见过不少同学上来就用实测数据结果滤波参数没调对稳定图一团糟最后怀疑代码有问题。先跑通仿真数据心里就有底了。1.3 SSI-COV和SSI-DATA怎么选代码包里两个版本都有刚开始我也纠结到底用哪个。SSI-COV协方差驱动先把输出数据的协方差矩阵算出来再构造Toeplitz矩阵做SVDSSI-DATA数据驱动则直接对Hankel矩阵做LQ分解本质上是在投影后对子空间进行估计。从使用体验上讲SSI-COV计算量小、内存占用低适合数据长度很长的情况这也是我默认推荐它做日常分析的原因。SSI-DATA理论上对噪声的鲁棒性更好数值稳定性更强但当数据量大时计算慢一些。这里给一个对比表对比项SSI-COVSSI-DATA核心操作协方差 Toeplitz SVDHankel LQ 投影 SVD计算速度快较慢内存占用低高数值稳定性良好更稳适合场景长数据、快速分析噪声大、数据较短实际项目中我的做法是先用SSI-COV跑一遍看稳定图效果。如果频率识别出来了但稳定轴不清晰再换SSI-DATA做交叉验证。两种方法的结果互相印证比单一方法更可信。2. 核心函数逐段拆解从原始数据到频率、阻尼比的全流程实现2.1 数据预处理去均值、去趋势、带通滤波预处理这一步常常被忽略但它对SSI结果的影响比算法本身还大。实测数据里传感器会有零漂结构本身可能受温度影响产生缓慢变形这些低频趋势在SSI眼里都是“伪模态”。如果不处理稳定图低频段会出现一串虚假极点特别迷惑人。我的PreProcess函数固定做三件事去均值、去线性趋势、带通滤波。function data PreProcess(data, fs, fLow, fHigh) % 去均值消除直流分量 data data - mean(data, 2); % 去线性趋势消除温度漂移等低频变化 data detrend(data, linear); % 带通滤波只保留关注频段 if ~isempty(fLow) ~isempty(fHigh) [b, a] butter(4, [fLow, fHigh] / (fs / 2), bandpass); data filtfilt(b, a, data); end end有几个细节要提一下。这里刻意用了filtfilt而不是filter因为filtfilt是零相位滤波不会产生相位偏移。普通filter会改变时间序列的相位特征对识别的阻尼比会造成明显偏差。数据矩阵统一用“通道数×采样点数”的格式这是SSI代码里最容易搞混的地方。滤波频段要根据实际关注范围定。比如桥梁环境振动通常关心0.2~10 Hz高频噪声较多就设到10 Hz。我的经验是带通范围宁可偏窄一点把无关频段滤干净稳定图的信噪比会有肉眼可见的提升但要保证目标模态频率都在带内。提示滤波阶数不要盲目太高。4阶巴特沃斯配合filtfilt已经足够阶数太高会引入边界振荡导致数据两端出现不可靠的响应段。2.2 协方差矩阵与Toeplitz矩阵的构造SSI-COV的起点是输出协方差。对于长度为N、通道数为l的响应数据y定义滞后量为k的协方差矩阵R(k) E[ y(tk) · y(t)^T ]这个矩阵反映了不同通道在不同时间滞后下的统计相关性。理论上结构模态信息就藏在这些协方差序列里。计算时可以写成两层循环但实际我会用矩阵乘法加速R zeros(l, l, i 1); for k 0:i R(:, :, k 1) (Y(:, k 1:end) * Y(:, 1:end - k)) / (size(Y, 2) - k); end这里的i是分块数它决定了后面Toeplitz矩阵的规模和可识别的最大系统阶数。接着把协方差序列排列成块Toeplitz矩阵T zeros(l * (i 1), l * (i 1)); for row 1:i 1 for col 1:i 1 T((row - 1) * l 1 : row * l, (col - 1) * l 1 : col * l) R(:, :, abs(row - col) 1); end endToeplitz矩阵的特点是每条对角线上的块都相同右下方不断重复R0、R1、R2……它把不同滞后量的协方差信息整合在一个矩阵里SSI后续的SVD分解实际上就是在对这个矩阵做低秩近似提取出系统可观测子空间。这里的分块数i不是随意取的。i太小Toeplitz矩阵高度不够高模态识别不出来i太大矩阵规模膨胀数值分解变慢且容易引入噪声。经验上i取目标最大阶数的一半到三分之二比较合适。比如预期系统在60阶以内i取3040即可。2.3 系统矩阵A、C的求解SVD与最小二乘有了块Toeplitz矩阵T之后随机子空间方法的核心数学操作就是奇异值分解。SVD能把T分解成U·S·V^T的形式其中奇异值大小代表各成分的“能量占比”。结构的真实模态对应较大的奇异值噪声对应较小的奇异值。function [A, C] SSI_COV(Y, i, n) % Y: l×N 响应数据 % i: 分块数 % n: 截断阶数一般为偶数 [l, N] size(Y); % 计算协方差 R zeros(l, l, i 1); for k 0:i R(:, :, k 1) (Y(:, k 1:end) * Y(:, 1:end - k)) / (N - k); end % 构造Toeplitz矩阵 T zeros(l * (i 1), l * (i 1)); for row 1:i 1 for col 1:i 1 T((row - 1) * l 1 : row * l, (col - 1) * l 1 : col * l) R(:, :, abs(row - col) 1); end end % SVD分解 [U, S, ~] svd(T, econ); % 截断到n阶 U1 U(:, 1:n); S1 S(1:n, 1:n); % 可观性矩阵 Oi U1 * sqrt(S1); % 利用移位结构求A A pinv(Oi(1:end - l, :)) * Oi(l 1:end, :); C Oi(1:l, :); end这个函数做了三件事算协方差、组Toeplitz、SVD截断后求A和C。看代码时要抓住一个关键点——n是系统阶数不是模态数量。一个模态对应一对共轭复数极点占2阶。因此如果期望识别3个模态n至少要取6实际稳定图分析中会在一定范围内扫描n。SVD截断后得到可观性矩阵Oi它满足一个重要的移位性质去掉最上面l行得到的矩阵等于去掉最下面l行的矩阵左乘A。用最小二乘就能解出A。C则直接取自Oi的前l行因为它对应输出方程中的观测矩阵。注意n必须是偶数否则提取出的特征值无法成对出现后续物理模态解释会很别扭。稳定图循环里一般从2开始、步长2递增。2.4 从系统矩阵到模态参数状态矩阵A是离散形式的系统矩阵它和连续时间系统之间有明确的关系。设采样间隔为dt 1/fs对A做特征值分解得到离散特征值λ对应的连续时间特征值s通过对数映射得到s ln(λ) / dt这样得到的s通常是一个复数实部代表衰减率虚部代表振荡角频率。欠阻尼结构的物理模态对应共轭复数对所以只取虚部为正的特征值即可。function [fn, damp, phi] ExtractModes(A, C, fs) dt 1 / fs; [Psi, D] eig(A); lambda diag(D); % 离散特征值 - 连续时间特征值 s_full log(lambda) / dt; % 只保留虚部为正的共轭对 idx imag(s_full) 0; s s_full(idx); Psi Psi(:, idx); % 固有频率 fn imag(s) / (2 * pi); % 阻尼比复平面中模长归一化 damp -real(s) ./ abs(s) * 100; % 振型输出矩阵乘特征向量 phi C * Psi; end阻尼比公式要特别注意。有人直接用-real(s)/imag(s)在阻尼比很小的时候两者差别不大但在阻尼比较大或复频率虚部与模长差别明显时除以abs(s)才更符合复平面几何关系。因为连续时间特征值s的模长就是系统固有圆频率ωn而实的衰减率是ξ·ωn相除自然得到阻尼比ξ。振型提取相对直接φ C·ψ其中ψ是A的特征向量。每一列对应一个模态幅值和相位信息都在里面。SSI得到的振型是复振型实际工程中如果结构阻尼不大复振型各分量的相位差应接近0或180度可以用这个特性判断识别出的“模态”是否物理。3. 实测流程参数配置、稳定图生成与算例验证3.1 几个关键参数的配置心得用好SSI参数配置比写代码本身更容易翻车。我把自己常用的参数和调试逻辑整理成一个表格方便直接对照参数常用范围调试说明采样率 fs由采集系统决定需高于目标最高频率的5~10倍以上滤波频带[fLow, fHigh]覆盖关注频段即可宁窄勿宽分块数 i30~60约为目标最大系统阶数的一半扫描阶数 n2:2:80从低到高遍历观察极点稳定性频率判稳阈值0.5%~1%相邻阶数极点频率相对偏差阻尼判稳阈值5%~10%阻尼比收敛较慢阈值可放宽MAC判稳阈值0.95以上振型一致性的硬指标分块数i是很多人最先踩坑的地方。想象一下Toeplitz矩阵的每一行都代表一个“观测窗口”窗口越多能容纳的系统阶数越高但窗口太多又会把噪声的统计特征也装进来。我一般在i40起步如果识别高频模态效果差逐步加大到60~80但很少超过100因为计算量会急剧上升。扫描阶数n的上限不要超过(i1)*l这是Toeplitz矩阵的实际行数。超过这个值SVD截断就失去了意义。通常nmax取60左右已经足够覆盖大多数土木结构的低阶模态。3.2 稳定图是怎么画出来的稳定图的核心思想很简单如果结构有真实模态那么无论系统阶数n取多少这个模态对应的极点在频率轴上应该差不多待在同一位置而噪声产生的虚假极点会随着n变化到处跑。把不同阶数下识别出的极点画在横轴频率、纵轴阶数的图上真实模态会形成一条竖直的“稳定轴”噪声极点则散落各处。StableDiagram.m的流程是for n nmin:2:nmax [A, C] SSI_COV(Y, i, n); [fn, damp, phi] ExtractModes(A, C, fs); % 与上一阶的极点做匹配 for k 1:length(fn) if ~isempty(prev_fr) [diff_min, j] min(abs(fn(k) - prev_fr)); if diff_min / prev_fr(j) 0.01 ... abs(damp(k) - prev_damp(j)) 10 ... MAC_calc(phi(:, k), prev_phi(:, j)) 0.95 % 稳定点红色标记 else % 新出现的点蓝色标记 end end end prev_fr fn; prev_damp damp; prev_phi phi; end这里有个提速技巧SSI稳定图本质上是不停地对同一个Toeplitz矩阵做不同截断阶数的SVD。与其每阶都重新做一次SVD不如先对T做一次完整SVD然后在循环里只做矩阵截取和最小二乘求解。代码包里的StableDiagram.m已经按这个思路优化了跑80阶也不会卡很久。如果你打算自己改写记得用这个方式。判稳条件里频率阈值最敏感一般设1%以内阻尼比识别稳定性天生较差阈值放宽到10%没问题MAC则要求0.95以上保证振型连续。有些文献会把“阻尼稳定”作为单独标记我在代码里选择同时满足三个条件才算“稳定点”这样稳定轴更干净代价是会漏掉一些阻尼比识别不稳定的真实模态。3.3 典型算例三自由度模拟结构的识别结果代码包自带的demo_data.mat是一个串联三自由度弹簧-质量系统的仿真数据。对它施加高斯白噪声激励采样频率100 Hz记录30分钟响应。理论上该系统前三阶固有频率为2.35 Hz、5.62 Hz、8.91 Hz阻尼比分别为1.8%、1.5%、1.2%。我用默认参数i40、n2:2:60跑了一遍稳定图上三条稳定轴非常清楚分别在2.3~2.4 Hz、5.6~5.7 Hz、8.9 Hz附近。取稳定轴上最密集的频率点作为最终识别结果模态理论频率/Hz识别频率/Hz频率误差/%理论阻尼比/%识别阻尼比/%第1阶2.352.3470.131.81.76第2阶5.625.6150.091.51.61第3阶8.918.9030.081.21.35频率识别精度非常高误差都在0.2%以内阻尼比误差在10%上下浮动。这个对比符合预期SSI对频率的估计十分可靠阻尼比受噪声和数据长度影响较大相同数据下结果会有一定波动。做实际项目时阻尼比建议多取几段数据平均或者结合对数衰减法做交叉验证不要只信单次识别的数值。4. 我踩过的坑常见问题与排查方法4.1 频率识别准、阻尼比却明显偏大这是SSI最常见的“症状”。我自己第一次跑实测数据时识别频率和有限元计算值对得很齐阻尼比却比设计值高了一倍一度怀疑代码有问题。排查后确认问题不在代码而在数据本身。环境激励的质量、数据长度、信噪比都会影响阻尼估计。激励不是理想白噪声、结构存在非线性、测试时间不够长这些因素会把额外的能量耗散“塞”进阻尼里导致识别值偏大。处理方法有几个增加数据长度。阻尼比估计对数据长度很敏感30分钟和5分钟的结果差距明显建议至少包含50个最低阶模态周期。分段平均。把数据切成多段分别识别再取平均能有效压制随机误差。放宽稳定图阻尼阈值。如果在判稳时要求阻尼比严格一致真实模态可能被当成不稳定点丢掉。与理论或有限元结果对比。如果阻尼比偏离预期太多不要急着调代码先怀疑数据是否满足白噪声激励假设。4.2 稳定图太乱没法选模态稳定图“毛刺”太多通常有三个原因滤波频带太宽、扫阶范围太高、判稳阈值太严。我调试时常采用“先粗后细”的策略。先把频率轴限定到关注范围比如只看0~20 Hz这样高频噪声极点会被排除在视野外。然后降低扫描阶数上限nmax从80降到40或30很多虚假极点会因为阶数不够而无法出现。最后如果还是乱就检查滤波——带通范围是否合适滤波阶数是否过高。还有一个容易被忽略的操作在生成稳定图前先看SSI_COV里SVD奇异值的曲线。奇异值会出现一个明显拐点拐点之前是信号主导拐点之后是噪声主导。这个拐点对应的大致阶数可以作为扫描上限的参考。4.3 运行报错与结果异常速查表整理了一些实际运行中容易遇到的问题按“现象→原因→处理办法”列成速查表现象可能原因处理方法Index exceeds array bounds数据长度N小于2i或通道数l设置错误检查输入数据矩阵尺寸确认是l×N格式特征值出现NaN或InfToeplitz矩阵奇异数据太短或通道冗余增大数据长度减少通道数或增加正则化提取的频率全是0滤波带通范围设置错误信号被滤除检查fLow/fHigh是否覆盖目标频率阻尼比偶尔出现负值对应的是计算极点而非物理模态这类极点直接舍弃不参与稳定图统计稳定图与上一阶完全一致阶数n超出Toeplitz矩阵实际秩SVD截断无意义降低nmax使其小于(i1)*l不同通道识别频率不一致通道间存在增益/相位不一致对传感器做标定检查同步采集是否正常内存不足i设置过大Toeplitz矩阵过大降低i或改用SSI_DATA方式分批处理注意SSI识别的结果里经常混入一些“数值模态”它们的阻尼比可能为负或异常大频率也不随阶数稳定。这些在稳定图上表现为孤立点不应该被选为最终模态。负阻尼单独出现时不要慌先检查它是不是在稳定轴上。最后再分享一个我从实践里总结的经验。很多人拿到SSI代码后第一件事就是往自己的实测数据上套结果一顿操作猛如虎稳定图却惨不忍睹。我的习惯是每到一个新项目先在仿真数据上把参数标定清楚再放到真实数据上微调滤波频带和扫描阶数。仿真数据的好处是理论值已知参数设置是否合理一目了然。这样做一轮之后再去处理实测数据你会有一种“稳了”的踏实感。本文还有配套的精品资源点击获取
