尧图网页设计 
简介基于SSI-COV算法的操作模态分析在MATLAB环境下形成了一套完整源码包面向结构工程与振动分析领域的科研人员和工程师用于识别结构在真实工况下的固有频率、模态形态与阻尼参数。压缩包共12个文件约2.39MB涵盖核心算法脚本、交互式示例、实测桥梁数据、说明文档及许可证等结构清晰便于按需查阅。目前已有52人学习适合结构工程背景的高年级学生、研究生及工程技术人员参考与二次开发。包内同时提供不依赖专业工具箱的独立算法版本与标准算法版本并附有稳定性分析绘图功能可辅助判断模态识别效果基于真实桥梁数据的交互式示例完整展示了从数据读取、协方差矩阵构建到模态参数提取的全流程有助于系统掌握SSI-COV方法的数学原理与MATLAB实现技巧。1. 操作模态分析场景下为什么选SSI-COV算法而不是峰值拾取在桥梁监测、风机塔筒和大型厂房改造这类工程里传感器测到的往往只有环境激励下的加速度时程激励力不可测。传统试验模态分析中依赖频响函数的方法彻底失效很多工程师的第一反应是将时程做FFT在幅值谱上找峰值。但频率密集、阻尼偏大或响应中存在拍频时峰值法给出的“模态”会在不同时间窗口漂移阻尼比更是没有可靠解释。SSI-COV算法的思路是绕过频谱直接从输出协方差中估计状态空间模型用SVD得到可观测矩阵再对系统矩阵做特征值分解。整个过程不选峰、不加窗、不依赖频响函数拟合对低频、密频和噪声都更稳健。下面用一套可运行的MATLAB代码串起从原始加速度时程到稳定图的完整链路。适合已经掌握振动力学基础、第一次自己写OMA程序的人也适合想验证第三方库识别结果的分析工程师。2. SSI-COV算法原理从随机响应到可观测矩阵2.1 时域算法与频域算法的边界操作模态分析的核心约束是只有输出y(t)输入激励无法记录。频域方法最容易实现的是以功率谱密度PSD为基础做奇异值分解FDD在SVD奇异值曲线上看到峰即模态频率。FDD的计算量小、参数少但有两个先天问题第一频率分辨率由FFT长度决定短数据频谱泄漏会掩盖密频模态第二阻尼是通过谱峰带宽拟合得到的对噪声和谱线间距极其敏感往往偏离真实值数倍。SSI-COV属于时域方法不把数据切成有限长度做频谱而是直接利用多个时刻响应之间的协方差因此不受FFT分辨率束缚。表1列出常用算法在OMA场景下的取舍。算法输入类型频率识别阻尼识别闭环适用主要代价FDD响应功率谱较好较差否需要平滑和峰值拾取SSI-DATA原始响应堆叠好好是QR分解开销大SSI-COV响应协方差好好是需要构造Toeplitz块矩阵ERA脉冲响应/互相关好较好是需要先得到IRF对于桥梁、风机这类以平稳环境激励为主的监测数据我一般优先选SSI-COV它对激励白噪声假设的要求比FDD宽松计算内存又比SSI-DATA小。MATLAB实现时只需要矩阵乘法和SVD不依赖System Identification Toolbox里的高阶封装函数。2.2 随机状态空间模型与协方差序列SSI-COV的起点是离散时间随机状态空间模型x(k1)A x(k)w(k)y(k)C x(k)v(k)其中x是状态向量y是ch通道的响应w、v为过程噪声和测量噪声统计上要求为白噪声且与当前状态无关。A是状态矩阵结构本身没有物理含义模态信息藏在A的特征值中这是后面需要辨识的目标。对输出协方差定义L(j)E[y(kj) y(k)^T]由于状态序列是平稳随机过程可以推出当j0时L(j)C A^(j-1) G。也就是说协方差序列中包含了与A、C一致的动力学信息。与数据驱动SSI-DATA不同SSI-COV不直接对原始响应堆叠矩阵做QR分解而是用不同滞后时间的协方差块构造一个分块Toeplitz矩阵T [L(0), L(-1), ..., L(-(i-1)); L(1), L(0), ..., L(-(i-2)); ... L(i-1), L(i-2), ..., L(0)]其中负滞后用L(-j)L(j)^T补齐。这个矩阵在理想条件下秩等于系统阶次n并可以分解为T O_i * Gamma_i其中O_i[C; C A; ...; C A^(i-1)]是扩展可观测矩阵Gamma_i是扩展可控性矩阵。只要能估计出O_i就能取出C和A。2.3 SVD分解是SSI-COV的识别核心对T做SVDT U S V^T。因为实际噪声导致矩阵满秩需要选择前r个主要奇异值并丢弃其余项。保留的奇异向量和奇异值重构出可观测矩阵的估计。MATLAB里的骨架代码如下[U,S,~] svd(T); O U(:,1:r) * sqrt(S(1:r,1:r)); C O(1:ch,:); A O(ch1:end,:) \ O(1:end-ch,:); [Psi,D] eig(A);从O_i的定义可知C就是前ch行。矩阵A利用可观测矩阵的移位结构把O_i去掉第一块得到O(ch1:end,:)把O_i去掉最后一块得到O(1:end-ch,:)两者满足O(ch1:end,:)O(1:end-ch,:)*A最小二乘即可解得A。得到A之后做特征值分解所有模态频率、阻尼比和振型都随之而来。整个识别链路里没有FFT、没有窄带滤波用户唯一需要调的关键参数是块行数i和截断阶次r这正是下一章用MATLAB落地时需要展开的部分。3. MATLAB实现数据预处理与块Toeplitz矩阵构造3.1 输入数据的最低要求和预处理先把数据统一成固定格式一个矩阵Y大小为通道数×采样点数行对应测点列对应时刻。采样率fs必须有准确值否则频率和阻尼无法换算成物理单位。对数据长度我的经验是最低模态周期的15~20倍以上比如最低关心频率0.5Hz至少需要30秒以上实际监测数据通常几分钟到几十分钟冗余度更高。通道数ch决定单次乘法的矩阵规模建议用6~16个测量通道过少会丢失空间振型过多则Toeplitz矩阵维数剧烈膨胀。预处理分三步走fs 256; % 采样率单位Hz Y detrend(Y, constant); % 先去掉均值 [bb, aa] butter(4, 40/(fs/2), low); % 4阶巴特沃斯低通 Y filtfilt(bb, aa, Y); % 零相位滤波避免相位偏移detrend第二个参数为constant时只去均值如果数据存在明显线性漂移可以改成默认的去线性趋势。低通滤波器上限设为关心频带上限的1~2倍比如此处40Hz。使用filtfilt而不是filter是因为零相位滤波不会让不同频率成分产生相对延迟代价是首尾约50个采样点失真这部分数据应在后续识别前裁掉。Y Y(:, 100:end-100); % 去掉filtfilt引入的首尾瞬态我一般按0.2~0.5秒计算裁掉的长度例如fs256时裁掉50~128点。若后续构建Toeplitz矩阵时块行数i较大数据起点再额外补裁i个点。预处理参数对识别结果的影响如下表参数推荐范围影响块行数 i20~50或i*ch达到模态数的3~5倍太小漏模态太大噪声极点增多滤波上限 fc1~2倍关心频带上限过高混入噪声过低损失模态数据长度最低模态周期的15~20倍以上不够则协方差估计方差大3.2 用MATLAB计算协方差序列与块Toeplitz矩阵协方差估计直接用矩阵乘法不需要循环通道function T build_ssicov_toeplitz(Y, i) % Y : 通道数×采样点数 % i : 输出矩阵的块行数同时也是单个Toeplitz块的尺寸 % T : 构造好的分块Toeplitz矩阵维数为(i*ch) x (i*ch) [ch, N] size(Y); % 1) 计算0~i阶滞后协方差 L zeros(ch, ch, i1); for k 0:i segLen N - k; L(:,:,k1) Y(:, k1:N) * Y(:, 1:segLen). / segLen; end % 2) 用滞后协方差铺满Toeplitz矩阵 T zeros(ch*i, ch*i); for r 1:i for c 1:i d r - c; if d 0 T((r-1)*ch1:r*ch, (c-1)*ch1:c*ch) L(:,:,d1); else T((r-1)*ch1:r*ch, (c-1)*ch1:c*ch) L(:,:,1-d).; end end end end说明L(:,:,k1)对应滞后k的协方差L(k)。segLenN-k是因为滞后k时只有N-k对乘积项可用。除以segLen是无偏估计实际用N也能跑差别在后几阶滞后上的估计噪声。铺矩阵时的关键是索引d0直接用第d1个协方差块d0时取滞后-d协方差再转置。这个转置一旦漏掉整个Toeplitz矩阵不再对称后续SVD得到的可观测矩阵会被破坏。块行数i是SSI-COV第一个必调参数。它既决定了能识别的最大阶次r_maxi*ch也决定了协方差最长滞后。i太小会漏掉低阶模态或大阻尼模态i太大会让高滞后阶段的协方差方差增大出现大量数值极点。我一般取i让i*ch约为预期主导模态数的3~5倍例如6通道、预期前10阶模态取i30~50。这个值不需要一次精确稳定图扫描会覆盖多种可能的r但i的上下限直接影响扫描范围。4. SVD截断、系统矩阵辨识与模态参数计算4.1 奇异值分解与截断阶次选择对上一章构造的T直接调用SVD[U, S, V] svd(T, econ); sig diag(S); % 方式1累积能量阈值 cumRatio cumsum(sig) / sum(sig); r find(cumRatio 0.99, 1, first); % 方式2按数量级落差选 % r length(sig) - sum(sig/max(sig)1e-6);SVD在实测数据中的专业含义是把协方差空间分成“可观测子空间”和“噪声子空间”。前r个奇异值对应真实的低维状态后面的奇异值主要来自测量噪声和数据长度有限。实际场景里奇异值曲线很少出现清晰断崖更多是连续下降所以我不建议只看累积能量阈值作为最终结论97%~99%都试一遍最后通过稳定图比较。上面的r只用做一次初筛稳定图会把r从2扫到i*ch-2所以这里截断是否精确并不致命。econ参数当T是方阵且满秩时结果与完全SVD一致只保留非零奇异值避免后续矩阵维数不必要的膨胀。U、V分别为左、右奇异向量矩阵。4.2 从可观测矩阵求系统矩阵A和C保留前r个奇异值重构可观测矩阵O U(:,1:r) * sqrt(S(1:r,1:r)); % O_i的精简估计 C O(1:ch, :); % 可观测矩阵第一块就是C Om O(ch1:end, :); % 去掉第一块对应C*A.. C*A^(i-1) Op O(1:end-ch, :); % 去掉最后一块对应C..C*A^(i-2) A Om \ Op; % 最小二乘解使得 Om*A约等于 Op这里Om是(i-1)*ch行Op也是解得A是r×r方阵。MATLAB的\对非方阵自动执行最小二乘效果等同pinv(Om)*Op但数值稳定性更好且不要求Om满秩。如果数据质量很差A会包含实部很大的特征值模态图里表现为频率极高或阻尼为负的散点后续稳定图会把这些点滤掉。还有一种常见写法是A Op \ Om取决于你定义可观测矩阵时A放在上移还是下移。本文按C; CA; CA^2...排列O(ch1:end,:)是O去掉第一块所以斜率关系是O(ch1:end,:) ≈ O(1:end-ch,:)*A不能写反。写反时A会变成它的逆对应模态频率全是负值这是最容易踩错的地方。4.3 离散特征值到连续时间模态参数得到的A是离散时间步上的状态矩阵特征值需要映射回连续时间dt 1/fs; [Psi, D] eig(A); lambda_d diag(D); % 从离散特征值到连续特征值 lambda_c log(lambda_d) / dt; % 频率Hz和阻尼比% fn abs(lambda_c) / (2*pi); zeta -real(lambda_c) ./ abs(lambda_c) * 100; % 振型物理坐标 C * 状态空间特征向量 Phi C * Psi;Psi是特征向量矩阵复数lambda_c是复数对应一对共轭特征值一个模态。fn取模是因为共轭对实部相同、虚部符号相反取模后得到同一个正频率。zeta是阻尼比公式中负号保证正阻尼输出正值如果出现负阻尼需要检查是不是特征值编号跨了共轭对。Phi的每一列是一个振型它包含幅值和相位实际输出时通常对每列做幅值归一化例如Phi(:,k)Phi(:,k)/max(abs(Phi(:,k)))。需要特别留意的单位陷阱fs是采样率dt就是秒。如果时间轴用的是毫秒或分钟频率会差10^3量级。日志里经常出现几百Hz甚至上千Hz的“模态”第一件事不是调参数而是检查fs与dt是否匹配。另外log(lambda_d)要求特征值不为零MATLAB中零特征值对应积分器或纯刚性模态应提前剔除。下面列出我在维护这段代码时遇到的高频误用误用现象fs单位错成kHz识别频率整体放大1000倍A的移位方向写反频率虚数、阻尼比负值未剔除零特征值出现0Hz附近的极低频散点在正式汇报前我建议用识别出的A、C重构协方差和实测协方差对比。取拟合值与实测值的相对误差作为整体模型质量指标误差在10%以内才算基本可信。这个验证步骤不额外要求工具箱直接对Phi、fn做协方差重建就能完成。5. 稳定图自动构建、参数调整与MATLAB工程化验证5.1 固定SVD、循环阶次的稳定图计算稳定图的原理很简单不给定单一阶次r而是让r从2逐步增加到i*ch-2对每个r重新提取一组模态。因为真实结构模态会在不同r里反复出现而噪声模态随机漂移所以将各阶结果画成“频率-阶次”散点图竖线聚集处就是物理模态。核心循环如下rList 2:2:ch*i-2; allModes []; for r rList Or U(:,1:r) * sqrt(S(1:r,1:r)); C_r Or(1:ch,:); Om_r Or(ch1:end,:); Op_r Or(1:end-ch,:); A_r Om_r \ Op_r; [Psi_r, D_r] eig(A_r); lam_d log(diag(D_r))/dt; for k 1:r if imag(lam_d(k)) 0 || abs(lam_d(k)) 1e-6 continue; end allModes [allModes; r, abs(lam_d(k))/(2*pi), ... -real(lam_d(k))/abs(lam_d(k))*100, k]; end end这里只在imag(lam)0时保留一个共轭半支避免同一个模态被重复计数。实测中同一个模态的微小频率漂移在±0.5%以内时视为稳定可以按频率排序后聚类。5.2 稳定阈值与自动聚类我常用的稳定判据如下表判据阈值频率偏差1%阻尼偏差5%MAC值0.98先以0.5Hz频率间隔为桶做直方图把峰值附近所有点合并为一个候选模态再对候选模态计算平均振型和某个参考振型做MAC验证。自动化流程可以用clusterdata(f_all, Cutoff, 0.01)跑层次聚类但注意输入频率幅值相差大应先把频率标准化到0~1。clusterdata的Cutoff对输出点密度敏感我通常先画直方图再按物理直觉微调聚类参数。5.3 三个值得检查的异常现象负阻尼连续出现最常见原因是数据段内激励不平稳带来状态估计偏差应对数据段加滑动窗口重算或换更长时间段。低频段出现固定宽峰往往不是真模态而是滤波前的低频漂移。确认预处理里detrend之后已经截掉首尾瞬态。MAC值长期低于0.9振型被测点布置遗漏应回到布置图检查是否有节点落在模态节点附近。建议在批处理代码注释里写清楚对所有候选模态先执行一次基于频率与MAC的稳定点筛选再用聚类替换人工读图。遇到恶劣数据时把i减半或把滤波器上限提高到3倍再跑一次往往比调任何阈值都管用。本文还有配套的精品资源点击获取 SSI-COV算法在操作模态分析中的原理与MATLAB实现
简介基于SSI-COV算法的操作模态分析在MATLAB环境下形成了一套完整源码包面向结构工程与振动分析领域的科研人员和工程师用于识别结构在真实工况下的固有频率、模态形态与阻尼参数。压缩包共12个文件约2.39MB涵盖核心算法脚本、交互式示例、实测桥梁数据、说明文档及许可证等结构清晰便于按需查阅。目前已有52人学习适合结构工程背景的高年级学生、研究生及工程技术人员参考与二次开发。包内同时提供不依赖专业工具箱的独立算法版本与标准算法版本并附有稳定性分析绘图功能可辅助判断模态识别效果基于真实桥梁数据的交互式示例完整展示了从数据读取、协方差矩阵构建到模态参数提取的全流程有助于系统掌握SSI-COV方法的数学原理与MATLAB实现技巧。1. 操作模态分析场景下为什么选SSI-COV算法而不是峰值拾取在桥梁监测、风机塔筒和大型厂房改造这类工程里传感器测到的往往只有环境激励下的加速度时程激励力不可测。传统试验模态分析中依赖频响函数的方法彻底失效很多工程师的第一反应是将时程做FFT在幅值谱上找峰值。但频率密集、阻尼偏大或响应中存在拍频时峰值法给出的“模态”会在不同时间窗口漂移阻尼比更是没有可靠解释。SSI-COV算法的思路是绕过频谱直接从输出协方差中估计状态空间模型用SVD得到可观测矩阵再对系统矩阵做特征值分解。整个过程不选峰、不加窗、不依赖频响函数拟合对低频、密频和噪声都更稳健。下面用一套可运行的MATLAB代码串起从原始加速度时程到稳定图的完整链路。适合已经掌握振动力学基础、第一次自己写OMA程序的人也适合想验证第三方库识别结果的分析工程师。2. SSI-COV算法原理从随机响应到可观测矩阵2.1 时域算法与频域算法的边界操作模态分析的核心约束是只有输出y(t)输入激励无法记录。频域方法最容易实现的是以功率谱密度PSD为基础做奇异值分解FDD在SVD奇异值曲线上看到峰即模态频率。FDD的计算量小、参数少但有两个先天问题第一频率分辨率由FFT长度决定短数据频谱泄漏会掩盖密频模态第二阻尼是通过谱峰带宽拟合得到的对噪声和谱线间距极其敏感往往偏离真实值数倍。SSI-COV属于时域方法不把数据切成有限长度做频谱而是直接利用多个时刻响应之间的协方差因此不受FFT分辨率束缚。表1列出常用算法在OMA场景下的取舍。算法输入类型频率识别阻尼识别闭环适用主要代价FDD响应功率谱较好较差否需要平滑和峰值拾取SSI-DATA原始响应堆叠好好是QR分解开销大SSI-COV响应协方差好好是需要构造Toeplitz块矩阵ERA脉冲响应/互相关好较好是需要先得到IRF对于桥梁、风机这类以平稳环境激励为主的监测数据我一般优先选SSI-COV它对激励白噪声假设的要求比FDD宽松计算内存又比SSI-DATA小。MATLAB实现时只需要矩阵乘法和SVD不依赖System Identification Toolbox里的高阶封装函数。2.2 随机状态空间模型与协方差序列SSI-COV的起点是离散时间随机状态空间模型x(k1)A x(k)w(k)y(k)C x(k)v(k)其中x是状态向量y是ch通道的响应w、v为过程噪声和测量噪声统计上要求为白噪声且与当前状态无关。A是状态矩阵结构本身没有物理含义模态信息藏在A的特征值中这是后面需要辨识的目标。对输出协方差定义L(j)E[y(kj) y(k)^T]由于状态序列是平稳随机过程可以推出当j0时L(j)C A^(j-1) G。也就是说协方差序列中包含了与A、C一致的动力学信息。与数据驱动SSI-DATA不同SSI-COV不直接对原始响应堆叠矩阵做QR分解而是用不同滞后时间的协方差块构造一个分块Toeplitz矩阵T [L(0), L(-1), ..., L(-(i-1)); L(1), L(0), ..., L(-(i-2)); ... L(i-1), L(i-2), ..., L(0)]其中负滞后用L(-j)L(j)^T补齐。这个矩阵在理想条件下秩等于系统阶次n并可以分解为T O_i * Gamma_i其中O_i[C; C A; ...; C A^(i-1)]是扩展可观测矩阵Gamma_i是扩展可控性矩阵。只要能估计出O_i就能取出C和A。2.3 SVD分解是SSI-COV的识别核心对T做SVDT U S V^T。因为实际噪声导致矩阵满秩需要选择前r个主要奇异值并丢弃其余项。保留的奇异向量和奇异值重构出可观测矩阵的估计。MATLAB里的骨架代码如下[U,S,~] svd(T); O U(:,1:r) * sqrt(S(1:r,1:r)); C O(1:ch,:); A O(ch1:end,:) \ O(1:end-ch,:); [Psi,D] eig(A);从O_i的定义可知C就是前ch行。矩阵A利用可观测矩阵的移位结构把O_i去掉第一块得到O(ch1:end,:)把O_i去掉最后一块得到O(1:end-ch,:)两者满足O(ch1:end,:)O(1:end-ch,:)*A最小二乘即可解得A。得到A之后做特征值分解所有模态频率、阻尼比和振型都随之而来。整个识别链路里没有FFT、没有窄带滤波用户唯一需要调的关键参数是块行数i和截断阶次r这正是下一章用MATLAB落地时需要展开的部分。3. MATLAB实现数据预处理与块Toeplitz矩阵构造3.1 输入数据的最低要求和预处理先把数据统一成固定格式一个矩阵Y大小为通道数×采样点数行对应测点列对应时刻。采样率fs必须有准确值否则频率和阻尼无法换算成物理单位。对数据长度我的经验是最低模态周期的15~20倍以上比如最低关心频率0.5Hz至少需要30秒以上实际监测数据通常几分钟到几十分钟冗余度更高。通道数ch决定单次乘法的矩阵规模建议用6~16个测量通道过少会丢失空间振型过多则Toeplitz矩阵维数剧烈膨胀。预处理分三步走fs 256; % 采样率单位Hz Y detrend(Y, constant); % 先去掉均值 [bb, aa] butter(4, 40/(fs/2), low); % 4阶巴特沃斯低通 Y filtfilt(bb, aa, Y); % 零相位滤波避免相位偏移detrend第二个参数为constant时只去均值如果数据存在明显线性漂移可以改成默认的去线性趋势。低通滤波器上限设为关心频带上限的1~2倍比如此处40Hz。使用filtfilt而不是filter是因为零相位滤波不会让不同频率成分产生相对延迟代价是首尾约50个采样点失真这部分数据应在后续识别前裁掉。Y Y(:, 100:end-100); % 去掉filtfilt引入的首尾瞬态我一般按0.2~0.5秒计算裁掉的长度例如fs256时裁掉50~128点。若后续构建Toeplitz矩阵时块行数i较大数据起点再额外补裁i个点。预处理参数对识别结果的影响如下表参数推荐范围影响块行数 i20~50或i*ch达到模态数的3~5倍太小漏模态太大噪声极点增多滤波上限 fc1~2倍关心频带上限过高混入噪声过低损失模态数据长度最低模态周期的15~20倍以上不够则协方差估计方差大3.2 用MATLAB计算协方差序列与块Toeplitz矩阵协方差估计直接用矩阵乘法不需要循环通道function T build_ssicov_toeplitz(Y, i) % Y : 通道数×采样点数 % i : 输出矩阵的块行数同时也是单个Toeplitz块的尺寸 % T : 构造好的分块Toeplitz矩阵维数为(i*ch) x (i*ch) [ch, N] size(Y); % 1) 计算0~i阶滞后协方差 L zeros(ch, ch, i1); for k 0:i segLen N - k; L(:,:,k1) Y(:, k1:N) * Y(:, 1:segLen). / segLen; end % 2) 用滞后协方差铺满Toeplitz矩阵 T zeros(ch*i, ch*i); for r 1:i for c 1:i d r - c; if d 0 T((r-1)*ch1:r*ch, (c-1)*ch1:c*ch) L(:,:,d1); else T((r-1)*ch1:r*ch, (c-1)*ch1:c*ch) L(:,:,1-d).; end end end end说明L(:,:,k1)对应滞后k的协方差L(k)。segLenN-k是因为滞后k时只有N-k对乘积项可用。除以segLen是无偏估计实际用N也能跑差别在后几阶滞后上的估计噪声。铺矩阵时的关键是索引d0直接用第d1个协方差块d0时取滞后-d协方差再转置。这个转置一旦漏掉整个Toeplitz矩阵不再对称后续SVD得到的可观测矩阵会被破坏。块行数i是SSI-COV第一个必调参数。它既决定了能识别的最大阶次r_maxi*ch也决定了协方差最长滞后。i太小会漏掉低阶模态或大阻尼模态i太大会让高滞后阶段的协方差方差增大出现大量数值极点。我一般取i让i*ch约为预期主导模态数的3~5倍例如6通道、预期前10阶模态取i30~50。这个值不需要一次精确稳定图扫描会覆盖多种可能的r但i的上下限直接影响扫描范围。4. SVD截断、系统矩阵辨识与模态参数计算4.1 奇异值分解与截断阶次选择对上一章构造的T直接调用SVD[U, S, V] svd(T, econ); sig diag(S); % 方式1累积能量阈值 cumRatio cumsum(sig) / sum(sig); r find(cumRatio 0.99, 1, first); % 方式2按数量级落差选 % r length(sig) - sum(sig/max(sig)1e-6);SVD在实测数据中的专业含义是把协方差空间分成“可观测子空间”和“噪声子空间”。前r个奇异值对应真实的低维状态后面的奇异值主要来自测量噪声和数据长度有限。实际场景里奇异值曲线很少出现清晰断崖更多是连续下降所以我不建议只看累积能量阈值作为最终结论97%~99%都试一遍最后通过稳定图比较。上面的r只用做一次初筛稳定图会把r从2扫到i*ch-2所以这里截断是否精确并不致命。econ参数当T是方阵且满秩时结果与完全SVD一致只保留非零奇异值避免后续矩阵维数不必要的膨胀。U、V分别为左、右奇异向量矩阵。4.2 从可观测矩阵求系统矩阵A和C保留前r个奇异值重构可观测矩阵O U(:,1:r) * sqrt(S(1:r,1:r)); % O_i的精简估计 C O(1:ch, :); % 可观测矩阵第一块就是C Om O(ch1:end, :); % 去掉第一块对应C*A.. C*A^(i-1) Op O(1:end-ch, :); % 去掉最后一块对应C..C*A^(i-2) A Om \ Op; % 最小二乘解使得 Om*A约等于 Op这里Om是(i-1)*ch行Op也是解得A是r×r方阵。MATLAB的\对非方阵自动执行最小二乘效果等同pinv(Om)*Op但数值稳定性更好且不要求Om满秩。如果数据质量很差A会包含实部很大的特征值模态图里表现为频率极高或阻尼为负的散点后续稳定图会把这些点滤掉。还有一种常见写法是A Op \ Om取决于你定义可观测矩阵时A放在上移还是下移。本文按C; CA; CA^2...排列O(ch1:end,:)是O去掉第一块所以斜率关系是O(ch1:end,:) ≈ O(1:end-ch,:)*A不能写反。写反时A会变成它的逆对应模态频率全是负值这是最容易踩错的地方。4.3 离散特征值到连续时间模态参数得到的A是离散时间步上的状态矩阵特征值需要映射回连续时间dt 1/fs; [Psi, D] eig(A); lambda_d diag(D); % 从离散特征值到连续特征值 lambda_c log(lambda_d) / dt; % 频率Hz和阻尼比% fn abs(lambda_c) / (2*pi); zeta -real(lambda_c) ./ abs(lambda_c) * 100; % 振型物理坐标 C * 状态空间特征向量 Phi C * Psi;Psi是特征向量矩阵复数lambda_c是复数对应一对共轭特征值一个模态。fn取模是因为共轭对实部相同、虚部符号相反取模后得到同一个正频率。zeta是阻尼比公式中负号保证正阻尼输出正值如果出现负阻尼需要检查是不是特征值编号跨了共轭对。Phi的每一列是一个振型它包含幅值和相位实际输出时通常对每列做幅值归一化例如Phi(:,k)Phi(:,k)/max(abs(Phi(:,k)))。需要特别留意的单位陷阱fs是采样率dt就是秒。如果时间轴用的是毫秒或分钟频率会差10^3量级。日志里经常出现几百Hz甚至上千Hz的“模态”第一件事不是调参数而是检查fs与dt是否匹配。另外log(lambda_d)要求特征值不为零MATLAB中零特征值对应积分器或纯刚性模态应提前剔除。下面列出我在维护这段代码时遇到的高频误用误用现象fs单位错成kHz识别频率整体放大1000倍A的移位方向写反频率虚数、阻尼比负值未剔除零特征值出现0Hz附近的极低频散点在正式汇报前我建议用识别出的A、C重构协方差和实测协方差对比。取拟合值与实测值的相对误差作为整体模型质量指标误差在10%以内才算基本可信。这个验证步骤不额外要求工具箱直接对Phi、fn做协方差重建就能完成。5. 稳定图自动构建、参数调整与MATLAB工程化验证5.1 固定SVD、循环阶次的稳定图计算稳定图的原理很简单不给定单一阶次r而是让r从2逐步增加到i*ch-2对每个r重新提取一组模态。因为真实结构模态会在不同r里反复出现而噪声模态随机漂移所以将各阶结果画成“频率-阶次”散点图竖线聚集处就是物理模态。核心循环如下rList 2:2:ch*i-2; allModes []; for r rList Or U(:,1:r) * sqrt(S(1:r,1:r)); C_r Or(1:ch,:); Om_r Or(ch1:end,:); Op_r Or(1:end-ch,:); A_r Om_r \ Op_r; [Psi_r, D_r] eig(A_r); lam_d log(diag(D_r))/dt; for k 1:r if imag(lam_d(k)) 0 || abs(lam_d(k)) 1e-6 continue; end allModes [allModes; r, abs(lam_d(k))/(2*pi), ... -real(lam_d(k))/abs(lam_d(k))*100, k]; end end这里只在imag(lam)0时保留一个共轭半支避免同一个模态被重复计数。实测中同一个模态的微小频率漂移在±0.5%以内时视为稳定可以按频率排序后聚类。5.2 稳定阈值与自动聚类我常用的稳定判据如下表判据阈值频率偏差1%阻尼偏差5%MAC值0.98先以0.5Hz频率间隔为桶做直方图把峰值附近所有点合并为一个候选模态再对候选模态计算平均振型和某个参考振型做MAC验证。自动化流程可以用clusterdata(f_all, Cutoff, 0.01)跑层次聚类但注意输入频率幅值相差大应先把频率标准化到0~1。clusterdata的Cutoff对输出点密度敏感我通常先画直方图再按物理直觉微调聚类参数。5.3 三个值得检查的异常现象负阻尼连续出现最常见原因是数据段内激励不平稳带来状态估计偏差应对数据段加滑动窗口重算或换更长时间段。低频段出现固定宽峰往往不是真模态而是滤波前的低频漂移。确认预处理里detrend之后已经截掉首尾瞬态。MAC值长期低于0.9振型被测点布置遗漏应回到布置图检查是否有节点落在模态节点附近。建议在批处理代码注释里写清楚对所有候选模态先执行一次基于频率与MAC的稳定点筛选再用聚类替换人工读图。遇到恶劣数据时把i减半或把滤波器上限提高到3倍再跑一次往往比调任何阈值都管用。本文还有配套的精品资源点击获取