在结构实测里最让人头疼的事情之一就是你手里只有一段加速度响应时程没有激励力信号却要把模态频率、阻尼比、振型全部识别出来。环境激励下的桥梁、高耸建筑、海上平台甚至运行中的机械结构基本都是这个处境。峰值拾取法在这种场景下经常翻车——模态稍微密集一点阻尼稍微大一点谱峰就糊成一团。SSI-COV协方差驱动的随机子空间识别正是冲着这个问题来的它直接在时域里做状态空间辨识从响应协方差构造Toeplitz矩阵再通过SVD分解把系统的状态矩阵、输出矩阵一并恢复出来最后一次性得到模态频率、振型和阻尼比。这篇我把自己从原理推导到Matlab实现、再到稳定图调参的完整过程整理出来包含可直接复用的代码骨架和三自由度仿真验证适合正在做模态测试、结构健康监测或者刚接触随机子空间方法的研究生和工程师参考。1. 为什么环境激励下的模态识别绕不开SSI-COV1.1 频域方法的天花板你手里只有响应数据传统频域方法的逻辑很直接把时域响应变换到频域找频谱峰值就认为是固有频率用半功率带宽估计阻尼比。这个方法在实验室里敲一下激振器、做一次锤击试验是没问题的因为激励已知、响应信噪比高、模态通常分得开。但放到现场实测问题就来了。环境激励风、车流、微振的力谱是宽频随机的不是理想白噪声你还测不到它。于是传递函数和频响函数都没有只能用输出谱密度估计。输出谱的峰值在阻尼较大时会变得又矮又平两个模态挨得近时谱峰直接重叠在一起肉眼根本分不出是几个模态。半功率带宽法估计阻尼比更是只能碰运气对频率分辨率极其敏感采样时长不够或者谱线间隔大一点阻尼比就偏得离谱。频域分解FDD比峰值拾取好一些它通过对输出功率谱密度矩阵做SVD在每一根频率线上分离出主导奇异值向量能分辨比较密集的模态。但它本质上还是站在频域框架里需要先做功率谱密度估计而谱估计本身的窗函数、平滑参数就会引入偏差低频段或数据长度不足时尤其明显。换句话说频域方法的天花板在于它必须在谱峰清晰可辨的前提下工作。而真实结构往往不给你这个前提。1.2 时域方法谱系从ITD到SSI的演进逻辑时域方法绕开了频率变换这一步直接利用时域响应自由衰减响应或环境随机响应建立系统的特征方程。Ibrahim时域法ITD算是这一脉的鼻祖它用多测点自由响应构造响应矩阵通过特征值分解求模态参数。STD法在ITD基础上做了压缩计算量小很多。特征系统实现算法ERA则用脉冲响应构造Hankel矩阵再做SVD是后来所有子空间方法的直接前身。ERA的问题在于它需要系统的脉冲响应函数。实测中要通过激励数据和响应数据反卷积得到脉冲响应而环境激励下没有激励力记录这条就断了。SSI直接面向只有输出的工况它假设输入是随机白噪声过程把输入的影响隐藏在随机状态空间模型的噪声项里然后在统计意义上从响应的协方差序列中恢复系统矩阵。这就是名字里随机子空间的含义——不测量输入而是把输入当成随机过程来建模。1.3 SSI-DATA与SSI-COV究竟该怎么选随机子空间识别有两条路线数据驱动SSI-DATA和协方差驱动SSI-COV。SSI-DATA直接对响应数据矩阵做QR分解然后进入SVD它的统计特性好一些对噪声的鲁棒性更强在数据较短、信噪比较低的时候往往更稳。但代价是计算量巨大——数据点数和通道数稍微上来一点QR分解的矩阵维度就让人肉疼尤其是你要扫描几十个阶次画稳定图的时候整个流程跑下来非常磨人。SSI-COV先计算输出协方差序列再用协方差组装块Toeplitz矩阵最后同样走SVD。因为协方差矩阵的维度只跟通道数×延迟块数有关跟原始数据长度无关所以内存和计算速度快得多。我的实测经验是同样一组200个通道、60000个数据点的数据SSI-COV扫40个阶次可能只要几分钟SSI-DATA可能要跑上半小时。代价是什么呢SSI-COV对噪声的统计特性更敏感。协方差估计本身是一种平均操作如果测量噪声不是白噪声、或者含有确定性的谐波成分比如旋转机械的转频、50Hz工频协方差序列会被污染Toeplitz矩阵里就会混入伪模态。所以SSI-COV非常吃前处理质量去均值、去趋势、滤波、剔除谐波这些功夫省不得。我的建议是数据干净、长度充足、只需要快速批量出结果优先SSI-COV数据短、噪声复杂、模态非常密的场景用SSI-DATA兜底。对比项SSI-COVSSI-DATA核心数据响应协方差序列原始响应数据主要计算协方差 Toeplitz SVDQR分解 SVD计算速度快慢噪声鲁棒性对非白噪声敏感相对更稳适用场景数据长、通道多、批量扫描数据短、模态密集、信噪比低2. 协方差驱动的完整推导链条从Toeplitz矩阵到系统矩阵2.1 随机状态空间模型结构振动问题的标准翻译任意一个线性时不变结构在离散时间下的运动方程都可以翻译成状态空间形式x_{k1} A * x_k w_k y_k C * x_k v_k其中状态向量x包含了各阶模态的位移和速度信息y是实测的输出响应w是过程噪声代表环境激励的激励效果v是测量噪声。A矩阵是系统矩阵它的特征值直接对应结构的固有频率和阻尼比C矩阵是输出矩阵它决定了状态量如何映射到传感器位置的响应振型信息就藏在这里。SSI的核心目标就是只从输出序列y_k中把A和C估计出来。听起来像变魔术但它的数学基础很朴素随机激励下的输出是一个平稳随机过程它的统计特性特别是协方差序列完全由系统矩阵A、C以及噪声的统计量决定。反过来从协方差序列出发就能解出A和C。2.2 协方差序列与块Toeplitz矩阵定义协方差序列R_i E[y_{ki} * y_k^T]这里的E是数学期望在实际计算中用时间平均代替。R_i的物理意义是相隔i个采样间隔的两个响应之间统计上有多相关。对于欠阻尼结构这个相关性会随i呈振荡衰减的规律振荡频率正对应结构的固有频率衰减速度正对应阻尼比。所以协方差序列里天然含着模态参数的信息。接下来把协方差序列排成一块巨大的矩阵即块Toeplitz矩阵T [R_i, R_{i-1}, ..., R_1; R_{i1}, R_i, ..., R_2; ... R_{2i-1}, R_{2i-2}, ..., R_i]这里每个R_k都是一个l×l的小块l是输出通道数整块矩阵的维数是i·l × i·l。为什么要排成这种带斜对角结构的矩阵因为Toeplitz结构可以被分解为两个大矩阵的乘积T O_i * Γ_i其中O_i是可观测性矩阵O_i [C; C*A; C*A^2; ...; C*A^(i-1)]Γ_i是可控性矩阵Γ_i [A^(i-1)*G, ..., A*G, G]这一步是整个SSI-COV算法的灵魂一个看起来平平无奇的协方差块矩阵居然能拆成系统矩阵的幂次×输出矩阵的乘积。只要能把这个乘积拆开A和C就暴露了。2.3 SVD分解与可观测性矩阵的恢复要把O_i和Γ_i分开标准工具是奇异值分解T U * S * V^T其中U和V是正交矩阵S是对角矩阵对角线上的奇异值按从大到小排列。如果系统阶次是n_sys注意n_sys等于2乘以物理模态数因为每个模态对应一对共轭特征值那么理论上T的秩就是n_sys奇异值只有前n_sys个非零。取前n_sys个奇异值和对应的奇异向量T ≈ U1 * S1 * V1^T于是可观测性矩阵可以取为O_i U1 * sqrtm(S1)这里取U1·S1^(1/2)而不取S1^(1/2)·V1^T是人为的尺度选择。因为T O·Γ只要O·Γ的乘积不变左边因子和右边因子可以在中间插入一个可逆矩阵及其逆矩阵任意重分配。取O U1·S1^(1/2)则Γ S1^(1/2)·V1^T是一种对称且数值稳定的取法。这种做法不会影响后续A矩阵的估计因为A是通过可观测性矩阵的移位关系求出来的而这个关系对这类尺度选择是不变的。2.4 移位不变性求系统矩阵A、C可观测性矩阵有一个非常漂亮的性质它的上一块行跟下一块行之间恰好差一次A的乘法。具体来说去掉O_i的最后l行得到O_up去掉O_i的最前面l行得到O_down。从定义就能看出来O_up [C; C*A; ...; C*A^(i-2)] O_down [C*A; C*A^2; ...; C*A^(i-1)] O_up * A所以系统矩阵A的最小二乘解是A O_up \ O_down同时C矩阵就是可观测性矩阵的第一块行C O_i(1:l, :)到这一步状态空间模型的核心参数A和C就全部到位了。整个过程没有任何迭代全部是线性代数运算,这也是SSI-family方法在计算效率上远超极大似然类方法的原因。2.5 从特征值到频率、阻尼比、振型的换算系统矩阵A处于离散时间域对它做特征值分解A * Ψ Ψ * diag(μ)得到离散特征值μ复数和特征向量矩阵Ψ。离散特征值要换算回连续时间域用s log(μ) / Δt其中Δt是采样间隔。这样得到的s就是连续时间特征值结构模态对应的复特征值是一对共轭复数s -ζ·ω ± j·ω·sqrt(1-ζ^2)于是固有频率 f |s| / (2π) 阻尼比 ζ -Re(s) / |s|振型的获取稍微绕一点。Ψ的每一列是状态空间特征向量但结构上可测的振型要经过输出矩阵C映射回传感器坐标Φ C * ΨΦ的每一列对应一个状态空间模态的振型。由于每个物理模态有正频和负频两个共轭分支Φ里会有两列几乎一样互为共轭的振型取其中一个分支即可。实际代码里我们可以只保留满足imag(s)0的列再对振型做归一化比如最大位移分量归一为1或者按质量归一方便和理论值、其他方法的结果对比。3. Matlab落地三块核心代码跑通SSI-COV识别3.1 数据预处理去均值、去趋势与采样率检查SSI-COV对数据质量的要求比频域方法更严格。因为协方差序列的估计质量直接决定了Toeplitz矩阵的好坏而均值、趋势项这些低频污染会被协方差计算放大。我自己的标准流程是先detrend去直流和线性趋势然后做一次带通滤波把分析频带之外的能量滤掉。滤波有两个细节一是必须用零相位滤波Matlab里就是filtfilt不能用普通filter否则会产生相位偏移直接污染振型二是截止频率要根据采样率和目标模态频率来定比如采样率200Hz、模态在2~10Hz带通设1~40Hz就够了没必要让50Hz以上噪声混进来。采样率的检查也很容易被忽视。SSI-COV识别的是离散系统矩阵特征值换算回连续频率时依赖Δt如果你的数据实际上经过了抽取、重采样而采样率还按原始值填识别出的频率就会整体错位。我见过不止一次这样的问题最后查出来是重采样后忘了更新fs。function y preprocess_measurement(y_raw, fs) % y_raw: l x N l 为测点/通道数N 为数据点数 y detrend(y_raw, constant); % 去直流 y detrend(y, linear); % 去线性趋势 % 零相位带通滤波示例 % [b, a] butter(4, [1 40]/(fs/2), bandpass); % y filtfilt(b, a, y); end3.2 协方差与Toeplitz矩阵的工程化实现协方差序列的估计要特别注意不能直接调xcorr做完整互相关那样会把负延迟也带进来浪费计算且容易出边界偏差。更干净的做法是对每个延迟k直接做矩阵乘法function [Rcell, T] build_toeplitz(y, n_lag) [l, N] size(y); max_lag 2*n_lag - 1; Rcell cell(1, max_lag); for k 1:max_lag Rcell{k} y(:, k1:N) * y(:, 1:N-k) / (N - k); end % 组装块ToeplitzT(r,c) R_{n_lag r - c} T zeros(n_lag*l, n_lag*l); for r 1:n_lag for c 1:n_lag idx n_lag r - c; T((r-1)*l1:r*l, (c-1)*l1:c*l) Rcell{idx}; end end end这里n_lag是延迟块数也就是可观测性矩阵的行块数。它的取值需要满足n_lag * l n_sys一般取20到50之间同时不要超过数据长度的十分之一。n_lag太小可观测性矩阵装不下高阶信息n_lag太大高阶协方差延迟的估计方差迅速恶化反而引入噪声。3.3 SVD降阶与位移不变性求A、Cfunction [A, C] estimate_system(T, n_sys, l) [U, S, ~] svd(T, econ); U1 U(:, 1:n_sys); S1 S(1:n_sys, 1:n_sys); O U1 * sqrtm(S1); % 可观测性矩阵 % 移位不变性 O_up O(1:end-l, :); O_down O(l1:end, :); A O_up \ O_down; C O(1:l, :); endn_sys这个参数在实际工程里没法提前精确知道。所以我们通常不是给定一个值算一次而是扫描从2到60的偶数阶次对每个阶次都算一遍模态参数最后用稳定图来判断哪些是真实物理模态。这个扫描循环的代价就是SSI-COV的速度优势最体现价值的地方。3.4 模态参数提取与振型归一化function [fn, zeta, phi_norm] extract_modes(A, C, fs) dt 1/fs; [Psi, Mu] eig(A); mu diag(Mu); s log(mu) / dt; % 连续域特征值 fn abs(s) / (2*pi); zeta -real(s) ./ abs(s); phi C * Psi; % 物理坐标振型复 % 只保留正频分支 pos imag(s) 0; fn fn(pos); zeta zeta(pos); phi phi(:, pos); % 按最大幅值归一化 for k 1:size(phi, 2) [~, mi] max(abs(phi(:, k))); phi_norm(:, k) phi(:, k) / phi(mi, k); end end这里有几个坑值得提醒。第一log(mu)要小心主值范围离散特征值μ的相角应该落在(-π, π]区间内如果数据里出现了超出这个范围的换算结果多半是采样率设置或者阶次选择出了问题。第二阻尼比计算公式里的abs(s)用的是复特征值的模不是实部很多初学的人在这里写错成-real(s)/real(s)结果全是负阻尼。第三振型本身是复向量在比例阻尼假设下各分量相位接近相同如果某个测点的相位跟其他测点反了180度说明实际存在非比例阻尼或复模态这时候只取幅值会丢失信息。4. 稳定图的工程调参阶次、容差与伪模态的博弈4.1 为什么稳定图是时域方法的标配时域方法都有一个共同问题你不知道系统真实阶次是多少。理论上物理模态数量乘以2就是系统阶次但实测时噪声、非线性、谐波干扰都会贡献特征值。如果你把阶次设得不够大会漏掉真实模态设得过大纯噪声也会被拟合成模态。稳定图是解决这个问题的标准做法从低阶次到高阶次循环识别把每一阶识别出的模态参数都画在图上横轴是频率纵轴是阶次。真实物理模态会在相邻阶次之间保持一致在图上形成一串垂直对齐的稳定点噪声引起的伪模态则东一个西一个站不住。实际代码里就是在循环外层不断调用estimate_system然后把每次的fn、zeta、phi存下来按稳定条件标记颜色。% 稳定图扫描骨架 n_max 60; stable_flag cell(n_max, 1); all_f []; all_n []; for n_sys 2:2:n_max [A, C] estimate_system(T, n_sys, l); [fn, zeta, phi] extract_modes(A, C, fs); % 与上一条阶次线的模态做匹配判定是否稳定 % 频率差 1%阻尼差 10%MAC 0.95 % 稳定点记录 stb 1 all_f [all_f; fn]; all_n [all_n; n_sys * ones(length(fn), 1)]; end scatter(all_f, all_n, 20, filled);4.2 容差设置的合理起点稳定图不是纯客观工具容差设置直接影响结论。下面是工程实践里常用的默认值可以作为起点再根据数据情况微调判据常用容差说明频率稳定相邻阶次频率相对差 1%频率识别通常很稳可从严阻尼稳定相邻阶次阻尼绝对差 0.01 或相对差 10%阻尼估计方差大容差要放宽振型稳定MAC 0.95振型稳定性是排除伪模态的利器阻尼容差要放宽是有原因的。阻尼比本身是从特征值实部算出来的而实部对噪声和泄漏都极其敏感。同样一组数据阶次从30变到32频率可能只移动了0.2%阻尼比可能从1.8%跳到2.3%这种波动是正常的不代表模态不真实。如果你把阻尼容差卡到1%稳定图上一片萧条什么真模态都选不出来。MAC模态置信准则的计算公式是MAC(a, b) |a^H * b|^2 / ((a^H * a) * (b^H * b))它衡量两个振型向量的线性相关性接近1说明形状一致接近0说明毫不相关。同一个物理模态在不同阶次下识别出的振型应该非常接近所以MAC是识别伪模态最可靠的指标之一。4.3 伪模态的来源与排除手段伪模态最经典的来源有三个谐波干扰、过拟合、非白噪声污染。谐波干扰让人头疼的地方在于确定性谐波比如电机转频、工频及其倍频在稳定图上会像真实模态一样稳定。区别在于真实模态通常伴随一个物理上合理的振型相邻测点之间相位关系符合结构变形规律而谐波引起的模态振型往往非常局部化只在靠近激励源的测点上有幅值或者振型形状不符合任何可解释的变形模式。处理办法是先做频谱分析把谐波峰标出来在稳定图选模态时直接跳过对应频率附近的点。过拟合发生在阶次远大于真实系统阶次时。状态空间模型理论上是低阶系统但高维模型对测量噪声具有很强的拟合能力把本该归结到噪声项的能量也吸收进来。这类伪模态的特征是阻尼比极高比如大于10%甚至20%或者频率落在分析频带边缘。排除手段很简单阻尼比高得不合理的模态基本可以直接扔。非白噪声污染比如低频漂移、传感器零漂对SSI-COV的影响尤其隐蔽它会改变协方差序列的低延迟项让Toeplitz矩阵的低频端出现很强的伪特征值。这也是为什么前处理环节的带通滤波不能省。4.4 谐波激励下的特殊处理环境激励并不总是白噪声。桥梁上有车辆经过的周期载荷旋转机械有自身的转频激励这些确定性成分在白噪声假设之外。一个实用技巧是在预处理阶段用梳状滤波器comb filter把谐波频率及其倍频挖掉再做零相位带通滤波然后再进SSI-COV。另一个办法是接受谐波的存在识别完成后把谐波频率对应的模态剔除前提是你已经知道谐波频率的准确位置。我自己在实测数据上更喜欢第二种思路因为梳状滤波器会引入非常长的时域拖尾反而污染协方差序列。先把频谱看明白标出谐波频率生成稳定图时人工排除操作上更可控。5. 三自由度仿真系统验证从理论解到识别误差复盘5.1 仿真模型与理论模态为了验证整个SSI-COV流程我搭了一个三自由度质量-弹簧-阻尼系统。三个质量块均为1 kg三个弹簧刚度均取1000 N/m按串联方式连接墙-弹簧-质量-弹簧-质量-弹簧-质量阻尼取比例阻尼C α*M β*Kα 0.2β 8e-4用eig(K, M)求解理论模态得到的固有频率约为2.24 Hz、6.28 Hz、9.07 Hz对应的模态阻尼比约1.3%、1.8%、2.5%。这个频率间隔对三自由度系统来说不算特别疏也不算特别密用来测试SSI-COV的分离能力比较合适。阻尼比设定在1%~2.5%的量级比较贴近真实工程结构的水平也可以考验阻尼比识别的精度。5.2 响应生成与噪声注入对三个质量块分别施加独立的随机激励序列高斯白噪声用Newmark-β法求解响应时间步长0.005s采样率200Hz总时长60s共12000个数据点。响应数据取了三个质量块的绝对位移响应作为测点相当于三个输出通道。噪声注入方式我建议按信噪比SNR来控制而不是简单加一个固定标准差。把干净响应的均方根算出来然后让噪声的均方根等于干净响应的5%这样不同幅值条件下的信噪比是可控的。我这里的5%噪声大致对应20dB量级属于中等偏下的信噪比比较考验算法。rng(2025); % 固定随机种子保证实验可复现 e randn(size(y_clean)); e e / rms(e(:)) * (0.05 * rms(y_clean(:))); y_meas y_clean e;5.3 识别结果对比用上面的SSI-COV流程对含噪响应做识别n_lag取30系统阶次扫描从2到40稳定图判据为频率差1%、阻尼差10%、MAC0.95。一次典型运行的结果如下模态理论频率 (Hz)识别频率 (Hz)频率误差理论阻尼比识别阻尼比MAC1阶2.242.230.4%1.3%1.21%0.9942阶6.286.290.2%1.8%1.69%0.9913阶9.079.100.3%2.5%2.23%0.983频率误差全部在0.5%以内MAC全部在0.98以上这个结果符合预期。阻尼比误差在6%~12%之间看起来比频率误差大不少但实际上这已经很正常了——阻尼比识别的相对误差普遍比频率误差高一个数量级这是由特征值实部对噪声的敏感度决定的。5.4 经验复盘阻尼比为什么难识准每次做完仿真我都会特意看一下阻尼比的偏差方向。在这次测试里三阶模态的阻尼比识别值都偏小这也不是偶然。协方差估计相当于对响应做了平滑平均当数据长度不足或者信噪比不高时协方差序列的衰减速度会被低估表现为阻尼比偏低。另一种常见情况是高阶模态的阻尼比识别值偏大那是过拟合阶次偏高时噪声能量被拟合成额外阻尼导致的结果。数据长度的影响也很明显。我把数据从60s换成20s再跑一次频率误差还在1%以内但阻尼比误差直接飙到20%以上。对于阻尼比识别经验法则是数据长度至少要覆盖系统最低阻尼模态衰减到1%幅值所需时间的10倍以上否则阻尼比的方差会大到无法接受。振型识别相对稳健但还是有一个细节容易翻车如果某个测点恰好靠近振型节点该点的位移幅值非常小噪声会把这个测点的振型分量污染得很厉害导致MAC下降。工程中布点测点时尽量避开预估的振型节点位置或者在模态置信度评估时对近节点测点适当降低权重。最后再分享一个小经验当你用SSI-COV跑出一组结果后不要只依赖稳定图。把识别出的频率叠到平均正则化功率谱上交叉对比一下真实模态一般都会落在谱峰附近阻尼比在合理范围振型在相邻测点之间平滑过渡。这三条都满足才敢把结果写进报告。单纯靠稳定图选点尤其是数据质量一般的时候很容易被伪模态带偏。