在工程测试里模态参数识别是个绕不开的活。你建了一个有限元模型算出了前几阶频率可实测结构到底是多少得靠锤击或环境激励数据来验证。传统的频域方法比如峰值拾取、频域分解在阻尼较大或者模态靠得近的时候往往分不清两个峰到底是一阶还是两阶。这时候我一般会换用SSI-COV随机子空间-协方差驱动法来处理。它的一个突出优点是直接从时域响应数据出发不用像频域法那样做快速傅里叶变换精度高而且能同时给出模态频率、振型和阻尼比。这篇博文就围绕多自由度系统模态参数识别把SSI-COV方法从数学原理到Matlab实现完整走一遍。我会把代码框架、关键参数怎么选、踩过的坑都交代清楚。适合正在做结构动力学课程设计、桥梁或机械结构的实测模态分析、以及想在新老算法之间做对比验证的研究生和工程师。只要你有两通道以上的加速度响应时程数据照着下面思路就能识别出比较靠谱的模态参数。1. 为什么在时域法里我会优先选SSI-COV结构模态参数识别方法大体分两派频域法和时域法。频域法如峰值拾取法、频域分解法、PolyMAX依赖测量频响函数。它的前提是激励已知、可以在结构上施加人工激励。可是对于大型结构例如人行桥、风力发电机塔筒锤击或激振器往往难以实施只能利用风、交通或地脉动作为自然激励。此时各输入点的激励力不可精确已知输入谱近似为宽带平坦的白噪声响应谱里包含结构全部模态信息。SSI-COV适合这类环境激励下的模态识别只用输出数据就能定阶并识别参数不需要激励力时程这正是随机子空间方法的“随机”二字的含义。与同类时域法如ITD、STD和特征系统实现算法ERA相比SSI-COV方法的稳健性更好。ERA虽然也能用自由响应数据做参数识别但如果处理的是环境激励下的平稳随机响应必须先做随机减量技术得到自由衰减响应中间多了一步预处理可能引入误差。SSI-COV则直接把原始响应数据生成Hankel矩阵再对协方差序列进行奇异值分解在噪声存在的情况下能够用统计手段压制干扰。这个思路其实和主成分分析类似把信号空间和噪声空间分开正交投影滤掉大部分非结构相关分量。另外一个选择SSI-COV的原因是编程逻辑非常清晰在Matlab里实现大概只需要六七个核心函数。从状态空间方程出发整个方法可以分为五步构造Hankel矩阵、估计输出协方差序列、形成Toeplitz矩阵、对Toeplitz矩阵做奇异值分解、从分解后的矩阵中提取系统状态矩阵。后续模态参数频率、阻尼比、振型只是对状态矩阵做特征值分解的结果。这种模块化的流程很适合教学演示也方便做代码调试——每一步我都能打印中间矩阵的尺寸和数值检查问题出在哪一步。近些年还出现了多级SSI、预知化SSI等变体实验中最常用的还是经典SSI-COV和SSI-DATA数据驱动。SSI-DATA对非线性振动和环境非平稳干扰的适应能力更强但计算量明显更大。如果只是做离线分析数据段长度几十万点SSI-COV的计算代价可以接受而且选对定阶方式后结果和SSI-DATA几乎一致因此很多开源工具箱如Matlab的MA工具箱、OpenModal默认提供SSI-COV选项不是没有道理。2. SSI-COV的核心数学推导从状态空间到Toeplitz矩阵要理解SSI-COV的处理流程得先回到动力学基本方程。对于一个n自由度线性时不变系统其运动方程离散化后可以写成状态空间形式x_{k1} A x_k w_ky_k C x_k v_k其中x_k是2n维状态向量包含位移和速度y_k是m维输出向量比如m个测点的加速度响应w_k和v_k分别是过程噪声和测量噪声假设为零均值白噪声。A是离散状态矩阵C是输出矩阵。从这一步可以看出识别模态参数的任务变成了在只知道输出y_k的情况下估计矩阵A和C。A的特征值和特征向量与系统的极点、模态振型之间有一一对应关系。关键在于不知道输入激励w_k系统是闭环的吗不是A和C仍然属于能观子空间的一部分可以用随机子空间理论证明输出协方差序列里已经包含了足够信息。设输出协方差矩阵序列为R_i E[y_{ki} y_k^T]。这个量不需要激励信息只取决于响应统计特性。把这些协方差序列排列成Toeplitz矩阵T_{1|i} [ R_i, R_{i-1}, ..., R_1; R_{i1}, R_i, ..., R_2; ...; R_{2i-1}, R_{2i-2}, ..., R_i ]这个Toeplitz矩阵的维度是“块行数乘输出通道数”乘以“块列数乘输出通道数”。例如有8个测点块行数和块列数各取20则T是160×160的矩阵。这时候对T做奇异值分解T U S V^T根据状态空间理论T的秩等于系统阶次2n。所以我们从S矩阵里找出明显非零的奇异值个数就能确定系统阶次。把S分解为S1和S2两部分S2对应噪声子空间可以进一步得到可观测性矩阵O_i U1·S1^(1/2)然后通过最小二乘拟合得到A和C。这一步是SSI-COV名字的来源——识别过程从协方差矩阵R_i出发而不是直接使用原始数据驱动。关于解析过程的更多说明实际操作中x_{k1}与y_k的互协方差矩阵G E[x_{k1} y_k^T]也会出现在中间推导里真正计算时并不需要显式估计x_k因为随机子空间算法会把状态序列消去只留下输出协方差。这一点让编程实现容易很多。理解整个算法我有一个比较好用的类比可以把Toeplitz矩阵想象成一张“多帧拼合的照片”。每一帧都是不同时间差下的输出相关值把这些帧叠成一个矩阵然后用奇异值分解做三维压缩。凡是结构模态引起的相关分量会在主奇异值里留下显著痕迹噪声引起的相关分量则分布均匀对应小奇异值。这就像是区分一张合照里的真实人物和随机噪点主成分是人物轮廓小成分是镜头噪点阈值一画就能分离。3. Matlab实现从响应数据到模态参数全流程3.1 生成测试数据用Newmark法模拟多自由度系统响应在拿到实测数据之前建议先用已知参数的多自由度系统生成仿真数据用来验证代码正确性。这是我会反复强调的做法。没有基准结果的代码只能叫“跑通”不能叫“验证”。以三自由度质量-弹簧-阻尼系统为例质量矩阵M、阻尼矩阵C和刚度矩阵K设定好后用Newmark-β法计算系统在白噪声激励下的位移、速度和加速度响应。这里的关键是激励力向量只作用在其中一个自由度上但响应在三个自由度上都能获取。然后用ode45或者Newmark法求解得到足够长的响应信号。采样频率建议设成系统最高频率的10倍以上比如系统最高固有频率为10 Hzfs设为100 Hz或200 Hz比较合适。数据点数至少取2^14 16384点越长的数据在统计意义上越能抑制协方差估计噪声。Matlab里生成响应的简略示例fs 100; % 采样频率 100 Hz N 30000; % 数据点数 M [2 0 0; 0 1 0; 0 0 0.5]; % 质量矩阵 K 1000 * [2 -1 0; -1 2 -1; 0 -1 1]; % 刚度矩阵 C 0.05 * M 0.02 * K; % 比例阻尼 F randn(3, N) * 10; % 白噪声激励 [y, v, x] newmark_beta(M, C, K, F, fs);注意这里加速度输出y是后续SSI-COV的输入。使用白噪声激励是为了和SSI-COV的模型假设匹配。如果你手头有实测数据就没有这一步但也要检查数据是否近似平稳、均值是否为零。线性趋势必须去除否则协方差估计会产生虚假的慢变分量。3.2 协方差序列与Toeplitz矩阵的代码实现SSI-COV算法第一步是估计输出协方差。由于只有有限长度数据协方差R_i可以直接用延时相关估计得到[m, N] size(y); % m 是通道数N是数据点数 y y - mean(y, 2); % 去除均值 maxlen 20; % 最大延时/块行数 RLag zeros(m, m, maxlen); for i 1:maxlen y1 y(:, 1:N-i); y2 y(:, i1:N); RLag(:, :, i) (y1 * y2) / (N - i); end这里RLag(:, :, i)就是R_i。为了提高协方差估计精度可以根据通道数m和块行数i调整归一化方法有些工具箱还会对数据先做互功率谱加权预白化处理。在Matlab实现中如果你发现协方差序列在高延时段数值翻转剧烈多半是数据长度不够或非平稳成分没有去干净。有了R_i之后构造Toeplitz矩阵iBlock 20; % 块行数 T zeros(m*iBlock, m*iBlock); for r 1:iBlock for c 1:iBlock lagIdx r - c iBlock; if lagIdx 1 lagIdx maxlen T((r-1)*m1:r*m, (c-1)*m1:c*m) RLag(:, :, lagIdx); end end end实际代码中可以直接调用Matlab自带的toeplitz函数把各块协方差排列成块Toeplitz矩阵。但要小心块和标量不完全一样必须自己组装块不能直接把标量向量丢给toeplitz。好多新手在这里翻车矩阵维度对不上后面全乱。3.3 奇异值分解与系统定阶对T做奇异值分解[U, S, V] svd(T); singular_vals diag(S);绘制对数奇异值曲线前2n个奇异值会明显大于后面部分。系统阶次2n怎么选一个经验法则是奇异值从某个位置开始骤降之后缓慢下降的部分属于噪声子空间。把奇异值从大到小排列选择保留的个数常见做法是观察相邻奇异值比值计算singular_vals(1:end-1) ./ singular_vals(2:end)比值突然变大的地方对应截断点。也可以通过预设最大阶次60然后配合稳定图来判断。取定阶次为order比如6阶对应三自由度系统的6个状态变量则order 6; U1 U(:, 1:order); S1 S(1:order, 1:order); V1 V(:, 1:order);接下来构造可观测性矩阵。有多种等价形式常用的是O1 U1 * sqrtm(S1)。注意S1必须是方阵且正定如果出现NaN很可能是奇异值小于零。由于数值舍入误差极小奇异值可能被算成微负值这种情况下取绝对值或直接忽略都行。在稳定性要求较高的场合还可以用balance函数对求解过程做数值平衡减少因矩阵病态导致的误差。3.4 提取系统状态矩阵与模态参数由可观测性矩阵O1计算状态矩阵A的方法是利用其位移结构记O1为2n×2n矩阵把它分成上块和下块则下块等于上块乘以AO_top O1(1:(order-m), :); O_bot O1(m1:order, :); A_est O_top \ O_bot;这里的A_est是离散状态矩阵。输出矩阵C_est直接取O1的前m行即可。对A_est做特征值分解[V_eig, D_eig] eig(A_est); lambda diag(D_eig);对于采样时间Δt 1/fs连续时间极点与离散特征值的关系为s log(lambda) / Δt。系统自然圆频率为|s|阻尼比为-cos(angle(s))。具体来说把每个复数极点分解成实部σ和虚部ω_df_n abs(s) / (2*pi)ξ -real(s) / abs(s)振型则从输出矩阵C_est乘以特征向量矩阵得到Phi C_est * V_eig;这里得到的Phi每一列是复振型幅值代表振型形状相位代表测点间相对相位。实测中如果发现某些测点的相位偏离0或π说明这些测点附近存在局部非线性或阻尼非比例效应这时候振型信息要格外小心解读。模态频率、阻尼比和振型提取完毕后可以设计一个比较函数把理论值仿真时已知和识别值画在同一张图上。如果频率残差小于1%、阻尼比残差小于0.5%说明代码实现完全正确。如果偏差大优先怀疑协方差延时个数不够或定阶错误。4. 稳定图判断虚假模态的利器不管用什么方法识别模态最容易让人头疼的问题是噪声产生的虚假模态。直接看奇异值截断只是一个粗糙的定阶方式更专业可靠的方法是画稳定图。稳定图的基本做法是逐次增加系统阶次2、4、6、8……对每个阶次都做一遍SSI-COV识别出一组模态频率、阻尼比、振型。然后把这些结果按频率值画在横轴上纵轴显示对应阶次。如果某一阶模态在连续多个阶次下频率和阻尼比变化很小就认为它是真实结构模态稳定图中这些点会连成竖直的线。而噪声模态会随机跳动无法稳定下来。Matlab实现简述orders 2:2:60; freq_candidates []; damp_candidates []; mode_candidates []; for ori orders [A_i, C_i] ssi_cov(y, ori, iBlock); [fn_i, xi_i, phi_i] get_modes(A_i, C_i, fs); freq_candidates [freq_candidates, fn_i]; damp_candidates [damp_candidates, xi_i]; end稳定图判断模态真实性的常用阈值以百分比形式如下表所示参数稳定判据相对变化频率小于 1%阻尼比小于 5%振型MAC值大于 0.95注意阻尼比的稳定阈值要放宽因为阻尼比识别本身方差较大尤其是低阻尼结构阻尼比的变异系数可能达到20%到30%。如果你看到某条稳定线上阻尼比从2%跳到2.2%不要急着判它为虚假模态需要结合振型和频率综合判断。画稳定图时还有一个细节横轴频率范围不要画到采样率一半那么大应聚焦到所关心的频带。比如系统固有频率集中在0-20 Hz横轴到25 Hz就够了否则高频噪声模态一大堆图看起来很乱。稳定图函数里可以加一个频带筛选fn_sel fn_i(fn_i 0.5 fn_i 25);按这个范围把结果填进图里画面会清晰得多。另外在剪切阻尼比较大的结构里稳定图上会看到明显的“频率漂移”。同一阶模态低阶次识别出的频率比高阶次稍低或稍高这时选哪个值我的做法是取稳定段中阶次居中区域的均值。因为阶次过低时子空间未被充分扩展阶次过高时数值噪声开始干扰重根分离中间段相对可靠。5. 常见问题与排查技巧实录5.1 阻尼比识别为负值怎么办SSI-COV识别出负阻尼比这个问题我遇到过很多次。多数情况下不是系统真的不稳定而是协方差矩阵估计噪声导致极点跑到右半平面。出现负阻尼时先检查以下几点数据是否去均值去除线性趋势协方差延时长度是否过短数据段是否包含非线性响应一种工程处理是如果只有个别阶次的阻尼比略负比如-0.3%直接舍弃该阶次结果因为它多半是虚假模态但如果连续多阶都出现负阻尼可能是测点布置导致振型在该频段不可观需要增加测点或改变传感器方向。还有一种办法是增大Toeplitz矩阵的块数。块数从20增加到30相当于利用更长时间的响应相关性有时可以明显改善阻尼比估计。代价是矩阵维度增大计算耗时增加但离线分析完全可接受。5.2 振型归一化与符号问题SSI-COV识别的振型是绝对尺度无关的因为输出协方差里丢失了激励幅度信息。比较振型时必须做归一化。最常见的是按最大幅值归一化把振型向量的绝对值最大元素调整为1。这样做的好处是方便和有限元分析结果比较。振型符号也可能出现整体翻转。某次识别出的第一阶振型是[1, -0.5, 0.3]理论值是[-1, 0.5, -0.3]实际上它们是一样的模态。比较时用模态置信因子MAC来判断振型相似程度不要直观比较每个元素的符号。Mac值计算公式是MAC |φ_a^H φ_b|^2 / ((φ_a^H φ_a)(φ_b^H φ_b))MAC大于0.9时认为两阶模态高度相关。Matlab里几行就能算出来收录到比较脚本里非常方便。5.3 计算量与内存优化SSI-COV的计算瓶颈主要在奇异值分解。如果测点数量较多比如32通道块行数40Toeplitz矩阵是1280×1280svd一次在普通笔记本上大约需要几秒稳定图画40个阶次就要一两分钟。这在可接受范围内。如果你想进一步提升速度可以改用稀疏SVD只计算前几十个奇异值[U, S, V] svds(T, 80); % 只计算前80个奇异值这一步对于峰值内存占用也有改善。另外在构造Toeplitz矩阵之前可以对输出数据做降采样。前提是目标模态频率远低于奈奎斯特频率比如结构前几阶频率在20 Hz以下采样率如果高达1000 Hz可以先把数据降到200 Hz既保留模态信息又大幅减少计算量。降采样前记得用抗混叠滤波器否则高频噪声混叠到低频区会把低阶模态的阻尼比估计搞乱。5.4 环境激励不满足白噪声假设怎么办SSI-COV的数学模型假定输入为白噪声。实际中环境激励往往不是纯白噪声比如人行荷载有步频峰值约2 Hz附近风荷载在低频段谱密度不平坦。如果激励谱存在尖峰这些尖峰会在稳定图上形成附加的“结构模态”容易被误判为真实模态。处理办法是对应已知激励特征频率的谱峰通常不与结构模态稳定线重合或者即便重合也无法通过频率和阻尼比双重稳定性检验。实际操作中只需要在稳定图判读时格外小心并结合有限元预分析结果筛选频段。如果实测数据非平稳明显比如存在车辆刹停等瞬态脉冲建议先对数据进行分段处理取平稳段再进行SSI-COV分析。平稳段的长度不能少于系统自由衰减最慢模态周期的10倍过低会严重影响阻尼比估计。6. 实操心得与扩展方向用SSI-COV做模态参数识别我个人的体会是这个方法的代码实现并不难难在如何解释结果。仿真数据很容易跑出漂亮结果但到实测数据时传感器噪声、局部非线性、温度变化等因素都会让结果变得扑朔迷离。建议在正式分析前先做一个简单的频域峰值拾取快速了解结构大致频率范围然后来设SSI-COV的频率筛选范围效率会高很多。在Matlab里实现时我习惯把整个流程封装成几个独立函数ssi_cov_core核心算法、plot_stability稳定图、extract_modal_params从状态矩阵提取模态参数、mac_compare振型对比。这样方便换一组数据直接复用脚本。对于更高要求的用户可以直接参考开源的免费工具包如MA工具箱或OpenModal以及Matlab社区里多个版本的SSI实现脚本但拿工具包之前最好先用仿真数据验证它们是否满足你的设定条件。对于模型修正或者健康监测方向的应用SSI-COV识别出的模态参数可以作为有限元模型修正的目标。识别出多阶模态后用优化算法调整有限元模型中某些物理参数使得理论频率和振型MAC值与SSI-COV结果一致。这个方法在城市桥梁健康监测中应用很广。模态频率随环境温度变化的问题可以通过长期监测建立回归补偿模型但阻尼比的变异性需要更多实测数据支撑。最后再分享一个我的小技巧如果你手头的数据只有一个测点可别直接用SSI-COV因为单通道数据无法建立完整的空间振型信息。至少需要两个测点才能区分同频的重根模态而且测点数量越多识别重根模态的能力越强。所以实测布点时在目标模态的振型节点附近尽量多布置几个测点这比单纯提高采样率更能提升识别效果。把这个方法吃透之后你会发现手里的Matlab代码不只是跑通而已它可以作为进一步研究自动化模态识别、深度学习辅助定阶、随机减量技术等多种方向的起点。先从三自由度系统开始验证再过渡到连续梁、板类结构逐步积累经验遇到问题的时候你也能更快定位到底出在哪一环。
SSI-COV随机子空间法：从数学原理到Matlab模态参数识别全流程
在工程测试里模态参数识别是个绕不开的活。你建了一个有限元模型算出了前几阶频率可实测结构到底是多少得靠锤击或环境激励数据来验证。传统的频域方法比如峰值拾取、频域分解在阻尼较大或者模态靠得近的时候往往分不清两个峰到底是一阶还是两阶。这时候我一般会换用SSI-COV随机子空间-协方差驱动法来处理。它的一个突出优点是直接从时域响应数据出发不用像频域法那样做快速傅里叶变换精度高而且能同时给出模态频率、振型和阻尼比。这篇博文就围绕多自由度系统模态参数识别把SSI-COV方法从数学原理到Matlab实现完整走一遍。我会把代码框架、关键参数怎么选、踩过的坑都交代清楚。适合正在做结构动力学课程设计、桥梁或机械结构的实测模态分析、以及想在新老算法之间做对比验证的研究生和工程师。只要你有两通道以上的加速度响应时程数据照着下面思路就能识别出比较靠谱的模态参数。1. 为什么在时域法里我会优先选SSI-COV结构模态参数识别方法大体分两派频域法和时域法。频域法如峰值拾取法、频域分解法、PolyMAX依赖测量频响函数。它的前提是激励已知、可以在结构上施加人工激励。可是对于大型结构例如人行桥、风力发电机塔筒锤击或激振器往往难以实施只能利用风、交通或地脉动作为自然激励。此时各输入点的激励力不可精确已知输入谱近似为宽带平坦的白噪声响应谱里包含结构全部模态信息。SSI-COV适合这类环境激励下的模态识别只用输出数据就能定阶并识别参数不需要激励力时程这正是随机子空间方法的“随机”二字的含义。与同类时域法如ITD、STD和特征系统实现算法ERA相比SSI-COV方法的稳健性更好。ERA虽然也能用自由响应数据做参数识别但如果处理的是环境激励下的平稳随机响应必须先做随机减量技术得到自由衰减响应中间多了一步预处理可能引入误差。SSI-COV则直接把原始响应数据生成Hankel矩阵再对协方差序列进行奇异值分解在噪声存在的情况下能够用统计手段压制干扰。这个思路其实和主成分分析类似把信号空间和噪声空间分开正交投影滤掉大部分非结构相关分量。另外一个选择SSI-COV的原因是编程逻辑非常清晰在Matlab里实现大概只需要六七个核心函数。从状态空间方程出发整个方法可以分为五步构造Hankel矩阵、估计输出协方差序列、形成Toeplitz矩阵、对Toeplitz矩阵做奇异值分解、从分解后的矩阵中提取系统状态矩阵。后续模态参数频率、阻尼比、振型只是对状态矩阵做特征值分解的结果。这种模块化的流程很适合教学演示也方便做代码调试——每一步我都能打印中间矩阵的尺寸和数值检查问题出在哪一步。近些年还出现了多级SSI、预知化SSI等变体实验中最常用的还是经典SSI-COV和SSI-DATA数据驱动。SSI-DATA对非线性振动和环境非平稳干扰的适应能力更强但计算量明显更大。如果只是做离线分析数据段长度几十万点SSI-COV的计算代价可以接受而且选对定阶方式后结果和SSI-DATA几乎一致因此很多开源工具箱如Matlab的MA工具箱、OpenModal默认提供SSI-COV选项不是没有道理。2. SSI-COV的核心数学推导从状态空间到Toeplitz矩阵要理解SSI-COV的处理流程得先回到动力学基本方程。对于一个n自由度线性时不变系统其运动方程离散化后可以写成状态空间形式x_{k1} A x_k w_ky_k C x_k v_k其中x_k是2n维状态向量包含位移和速度y_k是m维输出向量比如m个测点的加速度响应w_k和v_k分别是过程噪声和测量噪声假设为零均值白噪声。A是离散状态矩阵C是输出矩阵。从这一步可以看出识别模态参数的任务变成了在只知道输出y_k的情况下估计矩阵A和C。A的特征值和特征向量与系统的极点、模态振型之间有一一对应关系。关键在于不知道输入激励w_k系统是闭环的吗不是A和C仍然属于能观子空间的一部分可以用随机子空间理论证明输出协方差序列里已经包含了足够信息。设输出协方差矩阵序列为R_i E[y_{ki} y_k^T]。这个量不需要激励信息只取决于响应统计特性。把这些协方差序列排列成Toeplitz矩阵T_{1|i} [ R_i, R_{i-1}, ..., R_1; R_{i1}, R_i, ..., R_2; ...; R_{2i-1}, R_{2i-2}, ..., R_i ]这个Toeplitz矩阵的维度是“块行数乘输出通道数”乘以“块列数乘输出通道数”。例如有8个测点块行数和块列数各取20则T是160×160的矩阵。这时候对T做奇异值分解T U S V^T根据状态空间理论T的秩等于系统阶次2n。所以我们从S矩阵里找出明显非零的奇异值个数就能确定系统阶次。把S分解为S1和S2两部分S2对应噪声子空间可以进一步得到可观测性矩阵O_i U1·S1^(1/2)然后通过最小二乘拟合得到A和C。这一步是SSI-COV名字的来源——识别过程从协方差矩阵R_i出发而不是直接使用原始数据驱动。关于解析过程的更多说明实际操作中x_{k1}与y_k的互协方差矩阵G E[x_{k1} y_k^T]也会出现在中间推导里真正计算时并不需要显式估计x_k因为随机子空间算法会把状态序列消去只留下输出协方差。这一点让编程实现容易很多。理解整个算法我有一个比较好用的类比可以把Toeplitz矩阵想象成一张“多帧拼合的照片”。每一帧都是不同时间差下的输出相关值把这些帧叠成一个矩阵然后用奇异值分解做三维压缩。凡是结构模态引起的相关分量会在主奇异值里留下显著痕迹噪声引起的相关分量则分布均匀对应小奇异值。这就像是区分一张合照里的真实人物和随机噪点主成分是人物轮廓小成分是镜头噪点阈值一画就能分离。3. Matlab实现从响应数据到模态参数全流程3.1 生成测试数据用Newmark法模拟多自由度系统响应在拿到实测数据之前建议先用已知参数的多自由度系统生成仿真数据用来验证代码正确性。这是我会反复强调的做法。没有基准结果的代码只能叫“跑通”不能叫“验证”。以三自由度质量-弹簧-阻尼系统为例质量矩阵M、阻尼矩阵C和刚度矩阵K设定好后用Newmark-β法计算系统在白噪声激励下的位移、速度和加速度响应。这里的关键是激励力向量只作用在其中一个自由度上但响应在三个自由度上都能获取。然后用ode45或者Newmark法求解得到足够长的响应信号。采样频率建议设成系统最高频率的10倍以上比如系统最高固有频率为10 Hzfs设为100 Hz或200 Hz比较合适。数据点数至少取2^14 16384点越长的数据在统计意义上越能抑制协方差估计噪声。Matlab里生成响应的简略示例fs 100; % 采样频率 100 Hz N 30000; % 数据点数 M [2 0 0; 0 1 0; 0 0 0.5]; % 质量矩阵 K 1000 * [2 -1 0; -1 2 -1; 0 -1 1]; % 刚度矩阵 C 0.05 * M 0.02 * K; % 比例阻尼 F randn(3, N) * 10; % 白噪声激励 [y, v, x] newmark_beta(M, C, K, F, fs);注意这里加速度输出y是后续SSI-COV的输入。使用白噪声激励是为了和SSI-COV的模型假设匹配。如果你手头有实测数据就没有这一步但也要检查数据是否近似平稳、均值是否为零。线性趋势必须去除否则协方差估计会产生虚假的慢变分量。3.2 协方差序列与Toeplitz矩阵的代码实现SSI-COV算法第一步是估计输出协方差。由于只有有限长度数据协方差R_i可以直接用延时相关估计得到[m, N] size(y); % m 是通道数N是数据点数 y y - mean(y, 2); % 去除均值 maxlen 20; % 最大延时/块行数 RLag zeros(m, m, maxlen); for i 1:maxlen y1 y(:, 1:N-i); y2 y(:, i1:N); RLag(:, :, i) (y1 * y2) / (N - i); end这里RLag(:, :, i)就是R_i。为了提高协方差估计精度可以根据通道数m和块行数i调整归一化方法有些工具箱还会对数据先做互功率谱加权预白化处理。在Matlab实现中如果你发现协方差序列在高延时段数值翻转剧烈多半是数据长度不够或非平稳成分没有去干净。有了R_i之后构造Toeplitz矩阵iBlock 20; % 块行数 T zeros(m*iBlock, m*iBlock); for r 1:iBlock for c 1:iBlock lagIdx r - c iBlock; if lagIdx 1 lagIdx maxlen T((r-1)*m1:r*m, (c-1)*m1:c*m) RLag(:, :, lagIdx); end end end实际代码中可以直接调用Matlab自带的toeplitz函数把各块协方差排列成块Toeplitz矩阵。但要小心块和标量不完全一样必须自己组装块不能直接把标量向量丢给toeplitz。好多新手在这里翻车矩阵维度对不上后面全乱。3.3 奇异值分解与系统定阶对T做奇异值分解[U, S, V] svd(T); singular_vals diag(S);绘制对数奇异值曲线前2n个奇异值会明显大于后面部分。系统阶次2n怎么选一个经验法则是奇异值从某个位置开始骤降之后缓慢下降的部分属于噪声子空间。把奇异值从大到小排列选择保留的个数常见做法是观察相邻奇异值比值计算singular_vals(1:end-1) ./ singular_vals(2:end)比值突然变大的地方对应截断点。也可以通过预设最大阶次60然后配合稳定图来判断。取定阶次为order比如6阶对应三自由度系统的6个状态变量则order 6; U1 U(:, 1:order); S1 S(1:order, 1:order); V1 V(:, 1:order);接下来构造可观测性矩阵。有多种等价形式常用的是O1 U1 * sqrtm(S1)。注意S1必须是方阵且正定如果出现NaN很可能是奇异值小于零。由于数值舍入误差极小奇异值可能被算成微负值这种情况下取绝对值或直接忽略都行。在稳定性要求较高的场合还可以用balance函数对求解过程做数值平衡减少因矩阵病态导致的误差。3.4 提取系统状态矩阵与模态参数由可观测性矩阵O1计算状态矩阵A的方法是利用其位移结构记O1为2n×2n矩阵把它分成上块和下块则下块等于上块乘以AO_top O1(1:(order-m), :); O_bot O1(m1:order, :); A_est O_top \ O_bot;这里的A_est是离散状态矩阵。输出矩阵C_est直接取O1的前m行即可。对A_est做特征值分解[V_eig, D_eig] eig(A_est); lambda diag(D_eig);对于采样时间Δt 1/fs连续时间极点与离散特征值的关系为s log(lambda) / Δt。系统自然圆频率为|s|阻尼比为-cos(angle(s))。具体来说把每个复数极点分解成实部σ和虚部ω_df_n abs(s) / (2*pi)ξ -real(s) / abs(s)振型则从输出矩阵C_est乘以特征向量矩阵得到Phi C_est * V_eig;这里得到的Phi每一列是复振型幅值代表振型形状相位代表测点间相对相位。实测中如果发现某些测点的相位偏离0或π说明这些测点附近存在局部非线性或阻尼非比例效应这时候振型信息要格外小心解读。模态频率、阻尼比和振型提取完毕后可以设计一个比较函数把理论值仿真时已知和识别值画在同一张图上。如果频率残差小于1%、阻尼比残差小于0.5%说明代码实现完全正确。如果偏差大优先怀疑协方差延时个数不够或定阶错误。4. 稳定图判断虚假模态的利器不管用什么方法识别模态最容易让人头疼的问题是噪声产生的虚假模态。直接看奇异值截断只是一个粗糙的定阶方式更专业可靠的方法是画稳定图。稳定图的基本做法是逐次增加系统阶次2、4、6、8……对每个阶次都做一遍SSI-COV识别出一组模态频率、阻尼比、振型。然后把这些结果按频率值画在横轴上纵轴显示对应阶次。如果某一阶模态在连续多个阶次下频率和阻尼比变化很小就认为它是真实结构模态稳定图中这些点会连成竖直的线。而噪声模态会随机跳动无法稳定下来。Matlab实现简述orders 2:2:60; freq_candidates []; damp_candidates []; mode_candidates []; for ori orders [A_i, C_i] ssi_cov(y, ori, iBlock); [fn_i, xi_i, phi_i] get_modes(A_i, C_i, fs); freq_candidates [freq_candidates, fn_i]; damp_candidates [damp_candidates, xi_i]; end稳定图判断模态真实性的常用阈值以百分比形式如下表所示参数稳定判据相对变化频率小于 1%阻尼比小于 5%振型MAC值大于 0.95注意阻尼比的稳定阈值要放宽因为阻尼比识别本身方差较大尤其是低阻尼结构阻尼比的变异系数可能达到20%到30%。如果你看到某条稳定线上阻尼比从2%跳到2.2%不要急着判它为虚假模态需要结合振型和频率综合判断。画稳定图时还有一个细节横轴频率范围不要画到采样率一半那么大应聚焦到所关心的频带。比如系统固有频率集中在0-20 Hz横轴到25 Hz就够了否则高频噪声模态一大堆图看起来很乱。稳定图函数里可以加一个频带筛选fn_sel fn_i(fn_i 0.5 fn_i 25);按这个范围把结果填进图里画面会清晰得多。另外在剪切阻尼比较大的结构里稳定图上会看到明显的“频率漂移”。同一阶模态低阶次识别出的频率比高阶次稍低或稍高这时选哪个值我的做法是取稳定段中阶次居中区域的均值。因为阶次过低时子空间未被充分扩展阶次过高时数值噪声开始干扰重根分离中间段相对可靠。5. 常见问题与排查技巧实录5.1 阻尼比识别为负值怎么办SSI-COV识别出负阻尼比这个问题我遇到过很多次。多数情况下不是系统真的不稳定而是协方差矩阵估计噪声导致极点跑到右半平面。出现负阻尼时先检查以下几点数据是否去均值去除线性趋势协方差延时长度是否过短数据段是否包含非线性响应一种工程处理是如果只有个别阶次的阻尼比略负比如-0.3%直接舍弃该阶次结果因为它多半是虚假模态但如果连续多阶都出现负阻尼可能是测点布置导致振型在该频段不可观需要增加测点或改变传感器方向。还有一种办法是增大Toeplitz矩阵的块数。块数从20增加到30相当于利用更长时间的响应相关性有时可以明显改善阻尼比估计。代价是矩阵维度增大计算耗时增加但离线分析完全可接受。5.2 振型归一化与符号问题SSI-COV识别的振型是绝对尺度无关的因为输出协方差里丢失了激励幅度信息。比较振型时必须做归一化。最常见的是按最大幅值归一化把振型向量的绝对值最大元素调整为1。这样做的好处是方便和有限元分析结果比较。振型符号也可能出现整体翻转。某次识别出的第一阶振型是[1, -0.5, 0.3]理论值是[-1, 0.5, -0.3]实际上它们是一样的模态。比较时用模态置信因子MAC来判断振型相似程度不要直观比较每个元素的符号。Mac值计算公式是MAC |φ_a^H φ_b|^2 / ((φ_a^H φ_a)(φ_b^H φ_b))MAC大于0.9时认为两阶模态高度相关。Matlab里几行就能算出来收录到比较脚本里非常方便。5.3 计算量与内存优化SSI-COV的计算瓶颈主要在奇异值分解。如果测点数量较多比如32通道块行数40Toeplitz矩阵是1280×1280svd一次在普通笔记本上大约需要几秒稳定图画40个阶次就要一两分钟。这在可接受范围内。如果你想进一步提升速度可以改用稀疏SVD只计算前几十个奇异值[U, S, V] svds(T, 80); % 只计算前80个奇异值这一步对于峰值内存占用也有改善。另外在构造Toeplitz矩阵之前可以对输出数据做降采样。前提是目标模态频率远低于奈奎斯特频率比如结构前几阶频率在20 Hz以下采样率如果高达1000 Hz可以先把数据降到200 Hz既保留模态信息又大幅减少计算量。降采样前记得用抗混叠滤波器否则高频噪声混叠到低频区会把低阶模态的阻尼比估计搞乱。5.4 环境激励不满足白噪声假设怎么办SSI-COV的数学模型假定输入为白噪声。实际中环境激励往往不是纯白噪声比如人行荷载有步频峰值约2 Hz附近风荷载在低频段谱密度不平坦。如果激励谱存在尖峰这些尖峰会在稳定图上形成附加的“结构模态”容易被误判为真实模态。处理办法是对应已知激励特征频率的谱峰通常不与结构模态稳定线重合或者即便重合也无法通过频率和阻尼比双重稳定性检验。实际操作中只需要在稳定图判读时格外小心并结合有限元预分析结果筛选频段。如果实测数据非平稳明显比如存在车辆刹停等瞬态脉冲建议先对数据进行分段处理取平稳段再进行SSI-COV分析。平稳段的长度不能少于系统自由衰减最慢模态周期的10倍过低会严重影响阻尼比估计。6. 实操心得与扩展方向用SSI-COV做模态参数识别我个人的体会是这个方法的代码实现并不难难在如何解释结果。仿真数据很容易跑出漂亮结果但到实测数据时传感器噪声、局部非线性、温度变化等因素都会让结果变得扑朔迷离。建议在正式分析前先做一个简单的频域峰值拾取快速了解结构大致频率范围然后来设SSI-COV的频率筛选范围效率会高很多。在Matlab里实现时我习惯把整个流程封装成几个独立函数ssi_cov_core核心算法、plot_stability稳定图、extract_modal_params从状态矩阵提取模态参数、mac_compare振型对比。这样方便换一组数据直接复用脚本。对于更高要求的用户可以直接参考开源的免费工具包如MA工具箱或OpenModal以及Matlab社区里多个版本的SSI实现脚本但拿工具包之前最好先用仿真数据验证它们是否满足你的设定条件。对于模型修正或者健康监测方向的应用SSI-COV识别出的模态参数可以作为有限元模型修正的目标。识别出多阶模态后用优化算法调整有限元模型中某些物理参数使得理论频率和振型MAC值与SSI-COV结果一致。这个方法在城市桥梁健康监测中应用很广。模态频率随环境温度变化的问题可以通过长期监测建立回归补偿模型但阻尼比的变异性需要更多实测数据支撑。最后再分享一个我的小技巧如果你手头的数据只有一个测点可别直接用SSI-COV因为单通道数据无法建立完整的空间振型信息。至少需要两个测点才能区分同频的重根模态而且测点数量越多识别重根模态的能力越强。所以实测布点时在目标模态的振型节点附近尽量多布置几个测点这比单纯提高采样率更能提升识别效果。把这个方法吃透之后你会发现手里的Matlab代码不只是跑通而已它可以作为进一步研究自动化模态识别、深度学习辅助定阶、随机减量技术等多种方向的起点。先从三自由度系统开始验证再过渡到连续梁、板类结构逐步积累经验遇到问题的时候你也能更快定位到底出在哪一环。
