简介一套基于随机子空间法的MATLAB模态识别代码面向模式识别、特征提取与结构模态分析场景适合计算机、电子信息、数学等专业学生及工程技术人员用于课程设计、期末大作业和毕业设计算法通过随机选择原始数据的特征与样本点构建多个子空间模型以降低高维复杂度并提升噪声鲁棒性。压缩包共五十七个文件其中四十三个m文件构成主体算法与脚本十三个txt文件提供多年份案例数据另有一个GUI组件包用于界面支持整体仅一点零三兆字节。代码采用参数化编程参数可灵活修改注释详尽并附带可直接运行的案例数据工程结构包括随机子空间类、主控制器与图形界面模块兼容MATLAB 2014a、2019a、2024a等版本已有一百八十八人学习下载。通过阅读源码和运行数据可系统掌握随机子空间法的数据准备、子空间构建、模态识别与结果分析全流程并获得清晰的代码框架和调试思路对高维数据降维建模及工程应用具有直接参考价值。1. 随机子空间法模态识别从一段加速度时程里把模态参数“抽”出来结构健康监测和振动台试验里最常遇到的一个需求是手里只有传感器采集的加速度响应没荷激励信号怎么把这套结构的固有频率、阻尼比和振型识别出来环境激励下的模态识别有一整类方法随机子空间法Stochastic Subspace IdentificationSSI是其中工程适用性最好的一支它直接把输出信号写成状态空间形式通过矩阵投影和奇异值分解求系统矩阵的特征值与特征向量从而一次性得到频率、阻尼和振型。对 IT 背景的工程技术人员来说这套方法的核心其实不在力学而在矩阵分解和数值线性代数的功底。本文围绕标题中的随机子空间法模态识别给出一套能在 Matlab 里跑通的最小实现并把 Hankel 矩阵分块、投影矩阵、子空间维数、稳定图判据这些容易踩坑的环节逐个讲透。适合需要处理实测加速度数据、又不想依赖商业软件内置黑箱模块的工程师。2. 为什么选协方差驱动的随机子空间法与其他输出Only方法的边界2.1 频域法和时域法的本质差异环境激励模态识别基本分两大阵营。频域法以峰值拾取Peak Picking和频域分解FDD为代表思路直观把加速度时程做 FFT功率谱峰值对应的频率就是模态频率。这个做法在阻尼小、模态稀疏的结构上很好用但碰到密集模态或阻尼较大时就露馅了——两个频率接近的峰在频谱上重合峰值拾取会把两个模态当成一个。时域法里常见的有 Ibrahim 时域法ITD、特征系统实现算法ERA和随机子空间法。ITD 和 ERA 本质上都依赖自由衰减响应或脉冲响应而环境激励下的输出是持续的白噪声激励响应没有明显的自由衰减段需要先通过自然激励技术NExT把互相关函数当作自由衰减来用。随机子空间法走的是另一条路直接假设激励为白噪声把问题转化为输出序列的协方差矩阵结构识别。随机子空间法内部还有分支最常用的是协方差驱动 SSISSI-COV和数据驱动 SSISSI-DATA两者的差别在于把 Hankel 矩阵分成“过去”和“未来”两块后怎么构造投影——前者用输出协方差矩阵的 Toeplitz 结构后者用 LQ 分解后的投影矩阵。工程上 SSI-COV 更常见它不需要迭代计算量小而且对噪声的鲁棒性在绝大多数土木工程场景下够用。2.2 SSI-COV 的数学模型建立假设结构在白噪声激励下的离散状态空间方程为x(k1) A x(k) w(k) y(k) C x(k) v(k)其中 x 是 n 维状态向量y 是 m 维输出向量m 是传感器通道数w 是过程噪声v 是测量噪声。关键是这个模型里没有输入项激励的影响被并入了噪声。系统的模态参数就藏在 A 矩阵的特征值和 C 矩阵的列向量里。这就是随机子空间法最核心的巧妙之处噪声 v 的存在让直接用输出序列拟合 A 变得困难但如果我们构造输出的协方差矩阵白噪声假设下协方差只在零延迟处有值非零延迟处的协方差只由系统动态决定。于是把不同延迟的协方差块 R(i) E[y(ki) y^T(k)] 排成 Toeplitz 矩阵这个矩阵就包含了系统所有可观测信息。具体实现时用有限长度数据的估计值代替理论协方差R zeros(m, m, 2*i-1); for lag 1:2*i-1 Y1 Y(:, 1:end-lag); Y2 Y(:, 1lag:end); R(:, :, lag) (Y2 * Y1) / size(Y1, 2); end代码中i是 Hankel 矩阵的行块数后面会细讲如何取值m是通道数。这段代码把不同时间延迟的输出互协方差全部算出来作为后续 Toeplitz 矩阵的原材料。实际使用中数据长度至少要是2*i*m的几倍否则协方差估计的方差会很大。3. Matlab 实现随机子空间法模态识别的完整流程3.1 数据准备和 Hankel 矩阵构造第一步是把原始加速度时程组织成块 Hankel 矩阵。设总采样点数为 N通道数为 m选择行块数 i通常取最大分析阶数的 2 到 3 倍构造H [ y(0) y(1) ... y(N-2i) ] [ y(1) y(2) ... y(N-2i1) ] [ ... ... ... ... ] [ y(i-1) y(i) ... y(N-i) ]这个 2i 行、N-2i1 列的矩阵上半块是“过去输出”下半块是“未来输出”。用过去输出去预测未来输出预测的秩就决定了系统的阶数——这是随机子空间法能在噪声中识别模态的根本原因。在 Matlab 里构造这个矩阵的常见做法是用hankel函数function [T, Yp, Yf] build_hankel(Y, i) % Y: m x N 加速度时程矩阵 % i: 行块数 [m, N] size(Y); T zeros(2*i*m, N-2*i1); for k 1:i rows (k-1)*m1 : k*m; T(rows, :) Y(:, k : N-i-k1); T(rows i*m, :) Y(:, ik : N-k1); end Yp T(1:i*m, :); Yf T(i*m1:2*i*m, :); end这里Yp是过去输出块Yf是未来输出块。裁剪掉末尾的2*i-1个点是为了保证所有数据都能被充分利用。i的取值直接影响识别精度取太小高阶模态丢失取太大协方差估计的样本数减少。工程上 i 取 20 到 100 之间具体要看采样频率和关注的最高模态频率。3.2 协方差驱动的投影矩阵计算SSI-COV 的核心步骤是构造过去输出的协方差块和未来-过去的互协方差块然后组成块 Toeplitz 矩阵function [U, S, V] ssi_cov(Y, i, method) % method: cov 协方差驱动, data 数据驱动 [T, Yp, Yf] build_hankel(Y, i); [m, N] size(Y); R zeros(i*m, i*m); for k 1:i Rk (Yf((k-1)*m1:k*m, :) * Yp(1:end, :)) / size(Yp, 2); R((k-1)*m1:k*m, :) Rk; end T zeros(i*m, i*m); for k 1:i T(:, (k-1)*m1:k*m) R(:, (i-k)*m1:(i-k1)*m); end T (T T) / 2; % 对称化 [U, S, V] svd(T); end这段代码先计算了未来输出块中每一行块与全部过去输出之间的互协方差再排列成块 Toeplitz 矩阵。最后一步对称化是数值稳定性技巧——理论上 Toeplitz 矩阵在数据无限长时是对称的有限长度下会有微小不对称强制对称可以避免后续特征分解出现复数伪根。关于method参数实际应用中还有一种选择是数据驱动 SSISSI-DATA它的做法是把 Hankel 矩阵做 LQ 分解再对 L 矩阵的特定块做 SVD。SSI-DATA 对噪声的鲁棒性略好但计算量大约是 SSI-COV 的三到五倍。在传感器少于 16 通道的典型应用里SSI-COV 的速度优势不明显我会看数据质量决定——如果信噪比极差、伪模态比较多优先试 SSI-DATA常规场景直接用 SSI-COV。3.3 从 SVD 结果提取系统矩阵对 Toeplitz 矩阵做奇异值分解后系统的可观测性矩阵等于U * sqrt(S)的前 n 行其中 n 是系统阶数。然后利用可观测性矩阵的时间平移不变性求出系统矩阵 Afunction [A, C, Lambda, Phi] extract_modes(U, S, n, m) % n: 系统阶数通常取 2*模态数 Un U(:, 1:n); Sn S(1:n, 1:n); Obs Un * sqrt(Sn); % 可观测性矩阵下移 m 行与上移 m 行对比 Obs_up Obs(1:end-m, :); Obs_down Obs(m1:end, :); A pinv(Obs_up) * Obs_down; C Obs(1:m, :); % 特征值分解 [V, D] eig(A); Lambda diag(D); Phi C * V; end这里pinv是伪逆用最小二乘意义下求解。Obs_up对应可观测性矩阵去掉最后 m 行Obs_down对应去掉最前面 m 行两者的关系是Obs_down Obs_up * A这就是系统矩阵 A 的估计来源。由于 A 是在离散时间域里它的特征值 λ 需要换算成连续域的频率和阻尼比omega abs(log(Lambda)) / dt; xi -real(log(Lambda)) ./ abs(log(Lambda)); freq omega / (2*pi);这里要注意log取主值特征值的角度决定了识别的频率范围。如果采样率是 fs能识别的最大频率不超过 fs/2这是由采样定理决定的与随机子空间法本身无关。实际工程中采样率至少要是关注最高模态频率的 5 到 10 倍否则高频模态的辨识精度会很差。4. 关键参数怎么定稳定图、Hankel 行块数、系统阶数4.1 稳定图判据和实现随机子空间法最麻烦的问题是输入的系统阶数 n 不直观。取小了漏模态取大了产生大量伪模态。工程界的标准解法是稳定图——把 n 从 2 到 60 逐渐增大对每个 n 做一次识别把结果画在一张图上横轴是频率纵轴是阶数真实模态在不同阶数下频率和阻尼比都基本稳定伪模态则散射分布。稳定判据的常见定义是相邻两个阶数之间频率变化小于 1%阻尼比变化小于 5%MAC 值大于 0.95。MACModal Assurance Criterion是衡量两个振型向量相关性的指标function [freq, xi, phi, labels] stability_diagram(Y, i, fs, n_max) [T, ~, ~] build_hankel(Y, i); prev struct(); labels cell(n_max/2, 1); freq []; xi []; phi []; for n 2:2:n_max [U, S, ~] svd(T); [A, C, ~, Phi] extract_modes(U, S, n, size(Y,1)); Lambda eig(A); omega abs(log(Lambda)) / (1/fs); freq_n omega / (2*pi); xi_n -real(log(Lambda)) ./ abs(log(Lambda)); % 取共轭对的一半 [~, idx] sort(freq_n, ascend); freq_n freq_n(idx(1:end/2)); xi_n xi_n(idx(1:end/2)); Phi_n Phi(:, idx(1:end/2)); label_n strings(length(freq_n), 1); for j 1:length(freq_n) if ~isempty(prev.freq) [~, pidx] min(abs(prev.freq - freq_n(j))); if abs(prev.freq(pidx) - freq_n(j)) / freq_n(j) 0.01 ... abs(prev.xi(pidx) - xi_n(j)) 0.05 ... mac(prev.phi(:,pidx), Phi_n(:,j)) 0.95 label_n(j) s; else label_n(j) o; end end end prev.freq freq_n; prev.xi xi_n; prev.phi Phi_n; freq [freq; freq_n]; xi [xi; xi_n]; phi [phi; Phi_n]; labels{n/2} label_n; end end注意这里 SVD 是重新对 T 矩阵做的没有复用之前的分解结果。一个常见的性能优化是只做一次 SVD然后在不同 n 下取前 n 个奇异值和对应向量——因为 SVD 的前 n 列不随 n 变化。上面的代码为了清晰没做这个优化实际数据量大时可以按下述方式修改[U, S, ~] svd(T); % 只做一次 % 然后循环里直接取 Un U(:, 1:n)稳定图判据中的 MAC 值在低阶数时很容易大于 0.95因为状态向量维度低振型之间的分辨力不足所以看稳定图时不要只看连续两三阶的稳定点至少连续五个阶数都标注为“s”才值得标成真实模态。4.2 Hankel 行块数的实际取值策略Hankel 行块数 i 是随机子空间法里最容易被忽略的参数。它的作用有两面增大 i能覆盖更长的相关时间有利于识别低频模态但 i 太大Toeplitz 矩阵的维度变大而用于估计协方差的数据段变短估计方差升高。经验公式是i floor(N / (10 * m))其中 N 是采样总点数m 是通道数。取 10 倍冗余是安全值。另一个约束来自采样率i * dt至少要覆盖最低关注模态的一个周期以上否则协方差在相关时间窗内还没衰减到零识别结果会有偏。写成选择逻辑就是i_min ceil(min_freq_period / dt / 2); % 至少半个周期 i_max floor(N / (5 * m)); i clamp(round(i_opt), i_min, i_max);阻尼比越小的结构协方差衰减越慢需要的 i 越大。对于阻尼比 1% 的桥梁相关时间可能长达几十个周期这时候 i 可能要放到 100 以上。4.3 模态数目的自适应选择除了画稳定图人工判断还有几个自动化辅助手段。经典的 AIC赤池信息准则和 BIC贝叶斯信息准则在系统辨识里常用来选阶但在随机子空间法中直接使用的效果一般因为假设的噪声模型和实际环境激励不匹配。较可靠的做法是观察奇异值谱真正模态对应的奇异值在一个量级上噪声对应的奇异值呈缓慢下降的斜坡。把奇异值从大到小排列计算相邻比值找最大跳变位置function n_est estimate_order(S, tol) sv diag(S); sv_norm sv / sv(1); % 从后往前找第一个低于容差的点 below find(sv_norm tol, 1); if isempty(below) n_est length(sv); else n_est below - 1; end end这个方法的局限在于模态密集时奇异值没有明显跳变。工业界实际用得最多的是先跑一个较宽的阶数范围画稳定图再用聚类算法把稳定点聚成几个簇每个簇对应一个模态。Matlab 2023b 之后的系统辨识工具箱里ssest函数支持自动阶数选择但其核心依然是类似的平衡截断思路用不用工具箱不影响理解随机子空间法的行为。5. 三种典型工程场景的命令级速查5.1 将 CSV 加速度数据导入并执行随机子空间法模态识别实测数据最常见的落地格式是 CSV第一列时间戳后面每列一个测点。很多人在导入环节就出错常见是用readtable把时间戳读成了datetime类型导致后续矩阵运算直接报错。稳妥的做法是读取后强制转数值data readmatrix(accel_response.csv); t data(:, 1); Y data(:, 2:end); % 转成 m x N fs 1 / mean(diff(t)); Y detrend(Y, constant); % 去均值避免直流分量干扰 % 可选: 带通滤波到关注频段 f_low 0.5; f_high 20; [b, a] butter(4, [f_low f_high]/(fs/2), bandpass); Yf filtfilt(b, a, Y); i 40; % 根据上一章策略调整 [U, S, V] ssi_cov(Yf, i, cov); [freq, xi, phi, st] stability_diagram(Yf, i, fs, 60);这段代码里的detrend去均值是必须的因为环境激励下的加速度响应往往有微小的直流偏置如果不消除协方差估计会受到常数分量的污染。滤波则会引入边缘效应所以用了filtfilt做零相位滤波。注意滤波会改变数据的统计特性如果只是为了识别模态尽量不要滤波除非数据里有明显的工频干扰。5.2 与有限元仿真结果的交叉验证识别出模态参数后需要验证结果不是伪模态。常用做法是用模态置信准则 MAC 来比对识别振型和有限元计算振型phi_fe load(fe_mode_shapes.mat).phi_fe; % 两个振型矩阵按列匹配 MAC_matrix zeros(size(phi,2), size(phi_fe,2)); for k 1:size(phi,2) for l 1:size(phi_fe,2) num abs(phi(:,k) * phi_fe(:,l))^2; den (phi(:,k)*phi(:,k)) * (phi_fe(:,l)*phi_fe(:,l)); MAC_matrix(k,l) num / den; end end % MAC 值大于 0.8 的对角项对应匹配模态还有一个容易忽略的频率验证方法用半功率带宽法从功率谱独立估算阻尼比和随机子空间法的结果对比如果差超过 50%优先怀疑随机子空间法的参数设置。这个交叉验证能帮你发现是否把算法参数调得太激进。5.3 批处理多组数据的试探性搜索健康监测系统的数据是持续不断的每次都要人工看稳定图不现实。实际项目的常见做法是把稳定图变成批量处理维度files dir(data/*.csv); N length(files); freq_all cell(N, 1); for k 1:N data_k readmatrix(fullfile(files(k).folder, files(k).name)); Yk data_k(:, 2:end); freq_all{k} run_ssi(Yk, fs, 40); end % 聚合所有时段的频率识别结果用游程编码找重复出现的频率在这种批处理模式下稳定图的“稳定”不再要求连续阶数改进为连续时段上的重复性。需要补充的是环境激励的非平稳性会影响每个时间窗内协方差估计的一致性因此每个时间窗的数据长度至少要有 10 倍最低频率周期比如关注 0.5Hz 以上的模态窗长不应短于 100 秒含 2 倍余量否则每次识别的频率波动会大到难以收敛。6. 阻尼比识别偏大的原因与一种修正技巧随机子空间法识别阻尼比普遍存在偏大的问题尤其在低信噪比时。原因是测量噪声在协方差估计中贡献了一个正的偏置等价于给系统增加了人工阻尼。一个典型的修正方法是在识别前先做一次白化处理。白化滤波器只改变激励的频谱颜色不改变系统极点但能降低噪声对协方差矩阵的污染% 使用 AR 模型估计噪声颜色再对输出做白化 ar_order 20; a_coeff lpc(Y(1,:), ar_order); % 对第一个通道估计 AR 系数 Y_white filter(a_coeff, 1, Y);另一个实用技巧是用状态空间模型残差做修正。识别出 A 和 C 后重新生成预测输出残差里主要是测量噪声X_pred zeros(size(A,1), size(Y,2)); X_pred(:, 1) randn(size(A,1), 1); Y_pred zeros(size(Y)); for k 2:size(Y,2) X_pred(:, k) A * X_pred(:, k-1); Y_pred(:, k-1) C * X_pred(:, k-1); end residual Y - Y_pred; % 残差的标准差可以作为信噪比参考 noise_std std(residual, 0, 2);如果残差标准差和原始信号标准差差不多说明模型阶数取值过低识别结果基本不可信。这个检查比任何理论判据都快。还有一个容易被忽视的数据长度下限问题。很多工程上测了 10 分钟数据就认为足够但协方差驱动的随机子空间法里Toeplitz 矩阵的每个元素是从有限长数据估计出来的估计方差大致与有效样本数成反比。经验法则是有效样本数至少要是 Toeplitz 矩阵维度的 10 倍以上。设采样率 100Hz、通道 12 个、i 取 40Toeplitz 矩阵维度是 480那么最少需要 4800 个有效样本也就是 48 秒数据。考虑到协方差在延迟端的样本数更少120 秒以上比较安全。不满足这个条件时先降采样或减小 i 是最直接的做法。阻尼比修正这类数据处理技巧往往比换用更复杂的算法更能提升识别精度在实测数据处理时值得优先试一遍。本文还有配套的精品资源点击获取
随机子空间法模态识别：从加速度时程提取频率阻尼与振型的 Matlab 实现
简介一套基于随机子空间法的MATLAB模态识别代码面向模式识别、特征提取与结构模态分析场景适合计算机、电子信息、数学等专业学生及工程技术人员用于课程设计、期末大作业和毕业设计算法通过随机选择原始数据的特征与样本点构建多个子空间模型以降低高维复杂度并提升噪声鲁棒性。压缩包共五十七个文件其中四十三个m文件构成主体算法与脚本十三个txt文件提供多年份案例数据另有一个GUI组件包用于界面支持整体仅一点零三兆字节。代码采用参数化编程参数可灵活修改注释详尽并附带可直接运行的案例数据工程结构包括随机子空间类、主控制器与图形界面模块兼容MATLAB 2014a、2019a、2024a等版本已有一百八十八人学习下载。通过阅读源码和运行数据可系统掌握随机子空间法的数据准备、子空间构建、模态识别与结果分析全流程并获得清晰的代码框架和调试思路对高维数据降维建模及工程应用具有直接参考价值。1. 随机子空间法模态识别从一段加速度时程里把模态参数“抽”出来结构健康监测和振动台试验里最常遇到的一个需求是手里只有传感器采集的加速度响应没荷激励信号怎么把这套结构的固有频率、阻尼比和振型识别出来环境激励下的模态识别有一整类方法随机子空间法Stochastic Subspace IdentificationSSI是其中工程适用性最好的一支它直接把输出信号写成状态空间形式通过矩阵投影和奇异值分解求系统矩阵的特征值与特征向量从而一次性得到频率、阻尼和振型。对 IT 背景的工程技术人员来说这套方法的核心其实不在力学而在矩阵分解和数值线性代数的功底。本文围绕标题中的随机子空间法模态识别给出一套能在 Matlab 里跑通的最小实现并把 Hankel 矩阵分块、投影矩阵、子空间维数、稳定图判据这些容易踩坑的环节逐个讲透。适合需要处理实测加速度数据、又不想依赖商业软件内置黑箱模块的工程师。2. 为什么选协方差驱动的随机子空间法与其他输出Only方法的边界2.1 频域法和时域法的本质差异环境激励模态识别基本分两大阵营。频域法以峰值拾取Peak Picking和频域分解FDD为代表思路直观把加速度时程做 FFT功率谱峰值对应的频率就是模态频率。这个做法在阻尼小、模态稀疏的结构上很好用但碰到密集模态或阻尼较大时就露馅了——两个频率接近的峰在频谱上重合峰值拾取会把两个模态当成一个。时域法里常见的有 Ibrahim 时域法ITD、特征系统实现算法ERA和随机子空间法。ITD 和 ERA 本质上都依赖自由衰减响应或脉冲响应而环境激励下的输出是持续的白噪声激励响应没有明显的自由衰减段需要先通过自然激励技术NExT把互相关函数当作自由衰减来用。随机子空间法走的是另一条路直接假设激励为白噪声把问题转化为输出序列的协方差矩阵结构识别。随机子空间法内部还有分支最常用的是协方差驱动 SSISSI-COV和数据驱动 SSISSI-DATA两者的差别在于把 Hankel 矩阵分成“过去”和“未来”两块后怎么构造投影——前者用输出协方差矩阵的 Toeplitz 结构后者用 LQ 分解后的投影矩阵。工程上 SSI-COV 更常见它不需要迭代计算量小而且对噪声的鲁棒性在绝大多数土木工程场景下够用。2.2 SSI-COV 的数学模型建立假设结构在白噪声激励下的离散状态空间方程为x(k1) A x(k) w(k) y(k) C x(k) v(k)其中 x 是 n 维状态向量y 是 m 维输出向量m 是传感器通道数w 是过程噪声v 是测量噪声。关键是这个模型里没有输入项激励的影响被并入了噪声。系统的模态参数就藏在 A 矩阵的特征值和 C 矩阵的列向量里。这就是随机子空间法最核心的巧妙之处噪声 v 的存在让直接用输出序列拟合 A 变得困难但如果我们构造输出的协方差矩阵白噪声假设下协方差只在零延迟处有值非零延迟处的协方差只由系统动态决定。于是把不同延迟的协方差块 R(i) E[y(ki) y^T(k)] 排成 Toeplitz 矩阵这个矩阵就包含了系统所有可观测信息。具体实现时用有限长度数据的估计值代替理论协方差R zeros(m, m, 2*i-1); for lag 1:2*i-1 Y1 Y(:, 1:end-lag); Y2 Y(:, 1lag:end); R(:, :, lag) (Y2 * Y1) / size(Y1, 2); end代码中i是 Hankel 矩阵的行块数后面会细讲如何取值m是通道数。这段代码把不同时间延迟的输出互协方差全部算出来作为后续 Toeplitz 矩阵的原材料。实际使用中数据长度至少要是2*i*m的几倍否则协方差估计的方差会很大。3. Matlab 实现随机子空间法模态识别的完整流程3.1 数据准备和 Hankel 矩阵构造第一步是把原始加速度时程组织成块 Hankel 矩阵。设总采样点数为 N通道数为 m选择行块数 i通常取最大分析阶数的 2 到 3 倍构造H [ y(0) y(1) ... y(N-2i) ] [ y(1) y(2) ... y(N-2i1) ] [ ... ... ... ... ] [ y(i-1) y(i) ... y(N-i) ]这个 2i 行、N-2i1 列的矩阵上半块是“过去输出”下半块是“未来输出”。用过去输出去预测未来输出预测的秩就决定了系统的阶数——这是随机子空间法能在噪声中识别模态的根本原因。在 Matlab 里构造这个矩阵的常见做法是用hankel函数function [T, Yp, Yf] build_hankel(Y, i) % Y: m x N 加速度时程矩阵 % i: 行块数 [m, N] size(Y); T zeros(2*i*m, N-2*i1); for k 1:i rows (k-1)*m1 : k*m; T(rows, :) Y(:, k : N-i-k1); T(rows i*m, :) Y(:, ik : N-k1); end Yp T(1:i*m, :); Yf T(i*m1:2*i*m, :); end这里Yp是过去输出块Yf是未来输出块。裁剪掉末尾的2*i-1个点是为了保证所有数据都能被充分利用。i的取值直接影响识别精度取太小高阶模态丢失取太大协方差估计的样本数减少。工程上 i 取 20 到 100 之间具体要看采样频率和关注的最高模态频率。3.2 协方差驱动的投影矩阵计算SSI-COV 的核心步骤是构造过去输出的协方差块和未来-过去的互协方差块然后组成块 Toeplitz 矩阵function [U, S, V] ssi_cov(Y, i, method) % method: cov 协方差驱动, data 数据驱动 [T, Yp, Yf] build_hankel(Y, i); [m, N] size(Y); R zeros(i*m, i*m); for k 1:i Rk (Yf((k-1)*m1:k*m, :) * Yp(1:end, :)) / size(Yp, 2); R((k-1)*m1:k*m, :) Rk; end T zeros(i*m, i*m); for k 1:i T(:, (k-1)*m1:k*m) R(:, (i-k)*m1:(i-k1)*m); end T (T T) / 2; % 对称化 [U, S, V] svd(T); end这段代码先计算了未来输出块中每一行块与全部过去输出之间的互协方差再排列成块 Toeplitz 矩阵。最后一步对称化是数值稳定性技巧——理论上 Toeplitz 矩阵在数据无限长时是对称的有限长度下会有微小不对称强制对称可以避免后续特征分解出现复数伪根。关于method参数实际应用中还有一种选择是数据驱动 SSISSI-DATA它的做法是把 Hankel 矩阵做 LQ 分解再对 L 矩阵的特定块做 SVD。SSI-DATA 对噪声的鲁棒性略好但计算量大约是 SSI-COV 的三到五倍。在传感器少于 16 通道的典型应用里SSI-COV 的速度优势不明显我会看数据质量决定——如果信噪比极差、伪模态比较多优先试 SSI-DATA常规场景直接用 SSI-COV。3.3 从 SVD 结果提取系统矩阵对 Toeplitz 矩阵做奇异值分解后系统的可观测性矩阵等于U * sqrt(S)的前 n 行其中 n 是系统阶数。然后利用可观测性矩阵的时间平移不变性求出系统矩阵 Afunction [A, C, Lambda, Phi] extract_modes(U, S, n, m) % n: 系统阶数通常取 2*模态数 Un U(:, 1:n); Sn S(1:n, 1:n); Obs Un * sqrt(Sn); % 可观测性矩阵下移 m 行与上移 m 行对比 Obs_up Obs(1:end-m, :); Obs_down Obs(m1:end, :); A pinv(Obs_up) * Obs_down; C Obs(1:m, :); % 特征值分解 [V, D] eig(A); Lambda diag(D); Phi C * V; end这里pinv是伪逆用最小二乘意义下求解。Obs_up对应可观测性矩阵去掉最后 m 行Obs_down对应去掉最前面 m 行两者的关系是Obs_down Obs_up * A这就是系统矩阵 A 的估计来源。由于 A 是在离散时间域里它的特征值 λ 需要换算成连续域的频率和阻尼比omega abs(log(Lambda)) / dt; xi -real(log(Lambda)) ./ abs(log(Lambda)); freq omega / (2*pi);这里要注意log取主值特征值的角度决定了识别的频率范围。如果采样率是 fs能识别的最大频率不超过 fs/2这是由采样定理决定的与随机子空间法本身无关。实际工程中采样率至少要是关注最高模态频率的 5 到 10 倍否则高频模态的辨识精度会很差。4. 关键参数怎么定稳定图、Hankel 行块数、系统阶数4.1 稳定图判据和实现随机子空间法最麻烦的问题是输入的系统阶数 n 不直观。取小了漏模态取大了产生大量伪模态。工程界的标准解法是稳定图——把 n 从 2 到 60 逐渐增大对每个 n 做一次识别把结果画在一张图上横轴是频率纵轴是阶数真实模态在不同阶数下频率和阻尼比都基本稳定伪模态则散射分布。稳定判据的常见定义是相邻两个阶数之间频率变化小于 1%阻尼比变化小于 5%MAC 值大于 0.95。MACModal Assurance Criterion是衡量两个振型向量相关性的指标function [freq, xi, phi, labels] stability_diagram(Y, i, fs, n_max) [T, ~, ~] build_hankel(Y, i); prev struct(); labels cell(n_max/2, 1); freq []; xi []; phi []; for n 2:2:n_max [U, S, ~] svd(T); [A, C, ~, Phi] extract_modes(U, S, n, size(Y,1)); Lambda eig(A); omega abs(log(Lambda)) / (1/fs); freq_n omega / (2*pi); xi_n -real(log(Lambda)) ./ abs(log(Lambda)); % 取共轭对的一半 [~, idx] sort(freq_n, ascend); freq_n freq_n(idx(1:end/2)); xi_n xi_n(idx(1:end/2)); Phi_n Phi(:, idx(1:end/2)); label_n strings(length(freq_n), 1); for j 1:length(freq_n) if ~isempty(prev.freq) [~, pidx] min(abs(prev.freq - freq_n(j))); if abs(prev.freq(pidx) - freq_n(j)) / freq_n(j) 0.01 ... abs(prev.xi(pidx) - xi_n(j)) 0.05 ... mac(prev.phi(:,pidx), Phi_n(:,j)) 0.95 label_n(j) s; else label_n(j) o; end end end prev.freq freq_n; prev.xi xi_n; prev.phi Phi_n; freq [freq; freq_n]; xi [xi; xi_n]; phi [phi; Phi_n]; labels{n/2} label_n; end end注意这里 SVD 是重新对 T 矩阵做的没有复用之前的分解结果。一个常见的性能优化是只做一次 SVD然后在不同 n 下取前 n 个奇异值和对应向量——因为 SVD 的前 n 列不随 n 变化。上面的代码为了清晰没做这个优化实际数据量大时可以按下述方式修改[U, S, ~] svd(T); % 只做一次 % 然后循环里直接取 Un U(:, 1:n)稳定图判据中的 MAC 值在低阶数时很容易大于 0.95因为状态向量维度低振型之间的分辨力不足所以看稳定图时不要只看连续两三阶的稳定点至少连续五个阶数都标注为“s”才值得标成真实模态。4.2 Hankel 行块数的实际取值策略Hankel 行块数 i 是随机子空间法里最容易被忽略的参数。它的作用有两面增大 i能覆盖更长的相关时间有利于识别低频模态但 i 太大Toeplitz 矩阵的维度变大而用于估计协方差的数据段变短估计方差升高。经验公式是i floor(N / (10 * m))其中 N 是采样总点数m 是通道数。取 10 倍冗余是安全值。另一个约束来自采样率i * dt至少要覆盖最低关注模态的一个周期以上否则协方差在相关时间窗内还没衰减到零识别结果会有偏。写成选择逻辑就是i_min ceil(min_freq_period / dt / 2); % 至少半个周期 i_max floor(N / (5 * m)); i clamp(round(i_opt), i_min, i_max);阻尼比越小的结构协方差衰减越慢需要的 i 越大。对于阻尼比 1% 的桥梁相关时间可能长达几十个周期这时候 i 可能要放到 100 以上。4.3 模态数目的自适应选择除了画稳定图人工判断还有几个自动化辅助手段。经典的 AIC赤池信息准则和 BIC贝叶斯信息准则在系统辨识里常用来选阶但在随机子空间法中直接使用的效果一般因为假设的噪声模型和实际环境激励不匹配。较可靠的做法是观察奇异值谱真正模态对应的奇异值在一个量级上噪声对应的奇异值呈缓慢下降的斜坡。把奇异值从大到小排列计算相邻比值找最大跳变位置function n_est estimate_order(S, tol) sv diag(S); sv_norm sv / sv(1); % 从后往前找第一个低于容差的点 below find(sv_norm tol, 1); if isempty(below) n_est length(sv); else n_est below - 1; end end这个方法的局限在于模态密集时奇异值没有明显跳变。工业界实际用得最多的是先跑一个较宽的阶数范围画稳定图再用聚类算法把稳定点聚成几个簇每个簇对应一个模态。Matlab 2023b 之后的系统辨识工具箱里ssest函数支持自动阶数选择但其核心依然是类似的平衡截断思路用不用工具箱不影响理解随机子空间法的行为。5. 三种典型工程场景的命令级速查5.1 将 CSV 加速度数据导入并执行随机子空间法模态识别实测数据最常见的落地格式是 CSV第一列时间戳后面每列一个测点。很多人在导入环节就出错常见是用readtable把时间戳读成了datetime类型导致后续矩阵运算直接报错。稳妥的做法是读取后强制转数值data readmatrix(accel_response.csv); t data(:, 1); Y data(:, 2:end); % 转成 m x N fs 1 / mean(diff(t)); Y detrend(Y, constant); % 去均值避免直流分量干扰 % 可选: 带通滤波到关注频段 f_low 0.5; f_high 20; [b, a] butter(4, [f_low f_high]/(fs/2), bandpass); Yf filtfilt(b, a, Y); i 40; % 根据上一章策略调整 [U, S, V] ssi_cov(Yf, i, cov); [freq, xi, phi, st] stability_diagram(Yf, i, fs, 60);这段代码里的detrend去均值是必须的因为环境激励下的加速度响应往往有微小的直流偏置如果不消除协方差估计会受到常数分量的污染。滤波则会引入边缘效应所以用了filtfilt做零相位滤波。注意滤波会改变数据的统计特性如果只是为了识别模态尽量不要滤波除非数据里有明显的工频干扰。5.2 与有限元仿真结果的交叉验证识别出模态参数后需要验证结果不是伪模态。常用做法是用模态置信准则 MAC 来比对识别振型和有限元计算振型phi_fe load(fe_mode_shapes.mat).phi_fe; % 两个振型矩阵按列匹配 MAC_matrix zeros(size(phi,2), size(phi_fe,2)); for k 1:size(phi,2) for l 1:size(phi_fe,2) num abs(phi(:,k) * phi_fe(:,l))^2; den (phi(:,k)*phi(:,k)) * (phi_fe(:,l)*phi_fe(:,l)); MAC_matrix(k,l) num / den; end end % MAC 值大于 0.8 的对角项对应匹配模态还有一个容易忽略的频率验证方法用半功率带宽法从功率谱独立估算阻尼比和随机子空间法的结果对比如果差超过 50%优先怀疑随机子空间法的参数设置。这个交叉验证能帮你发现是否把算法参数调得太激进。5.3 批处理多组数据的试探性搜索健康监测系统的数据是持续不断的每次都要人工看稳定图不现实。实际项目的常见做法是把稳定图变成批量处理维度files dir(data/*.csv); N length(files); freq_all cell(N, 1); for k 1:N data_k readmatrix(fullfile(files(k).folder, files(k).name)); Yk data_k(:, 2:end); freq_all{k} run_ssi(Yk, fs, 40); end % 聚合所有时段的频率识别结果用游程编码找重复出现的频率在这种批处理模式下稳定图的“稳定”不再要求连续阶数改进为连续时段上的重复性。需要补充的是环境激励的非平稳性会影响每个时间窗内协方差估计的一致性因此每个时间窗的数据长度至少要有 10 倍最低频率周期比如关注 0.5Hz 以上的模态窗长不应短于 100 秒含 2 倍余量否则每次识别的频率波动会大到难以收敛。6. 阻尼比识别偏大的原因与一种修正技巧随机子空间法识别阻尼比普遍存在偏大的问题尤其在低信噪比时。原因是测量噪声在协方差估计中贡献了一个正的偏置等价于给系统增加了人工阻尼。一个典型的修正方法是在识别前先做一次白化处理。白化滤波器只改变激励的频谱颜色不改变系统极点但能降低噪声对协方差矩阵的污染% 使用 AR 模型估计噪声颜色再对输出做白化 ar_order 20; a_coeff lpc(Y(1,:), ar_order); % 对第一个通道估计 AR 系数 Y_white filter(a_coeff, 1, Y);另一个实用技巧是用状态空间模型残差做修正。识别出 A 和 C 后重新生成预测输出残差里主要是测量噪声X_pred zeros(size(A,1), size(Y,2)); X_pred(:, 1) randn(size(A,1), 1); Y_pred zeros(size(Y)); for k 2:size(Y,2) X_pred(:, k) A * X_pred(:, k-1); Y_pred(:, k-1) C * X_pred(:, k-1); end residual Y - Y_pred; % 残差的标准差可以作为信噪比参考 noise_std std(residual, 0, 2);如果残差标准差和原始信号标准差差不多说明模型阶数取值过低识别结果基本不可信。这个检查比任何理论判据都快。还有一个容易被忽视的数据长度下限问题。很多工程上测了 10 分钟数据就认为足够但协方差驱动的随机子空间法里Toeplitz 矩阵的每个元素是从有限长数据估计出来的估计方差大致与有效样本数成反比。经验法则是有效样本数至少要是 Toeplitz 矩阵维度的 10 倍以上。设采样率 100Hz、通道 12 个、i 取 40Toeplitz 矩阵维度是 480那么最少需要 4800 个有效样本也就是 48 秒数据。考虑到协方差在延迟端的样本数更少120 秒以上比较安全。不满足这个条件时先降采样或减小 i 是最直接的做法。阻尼比修正这类数据处理技巧往往比换用更复杂的算法更能提升识别精度在实测数据处理时值得优先试一遍。本文还有配套的精品资源点击获取
