简介面向结构动力学与健康监测研究者的数据驱动随机子空间识别DDRSI源码包聚焦模态自动识别可从含噪声的实测响应中直接提取频率、阻尼比和振型无需建立先验理论模型适用于桥梁、建筑、航空航天器等基础设施的振动监测与参数辨识。压缩包共24个文件包括13个MAT数据文件与11个M脚本数据文件涵盖多组仿真与实测信号如ASCE基准模型、连续6自由度结构等脚本则覆盖数据分块矩阵构建、稳定图绘制、模态参数自动筛选等关键环节便于对照学习和二次开发。已有1116人学习下载。资源包含完整可运行的主程序与配套函数并附有详细算例可帮助使用者系统掌握数据驱动随机子空间的建模步骤、去噪策略与虚假模态剔除方法为实际工程结构的状态评估和健康监测提供可复用的代码基础。1. 随机子空间识别不激励、不赋初值从实测响应里把模态拆出来一台风机机舱上的加速度计连续记录了几天环境振动数据没人敲它、也没开激振器却要报出结构前三阶频率、阻尼比和振型——这是运行模态分析里的典型任务。随机子空间识别SSI是这类任务里用得最广的算法族之一而数据驱动版本SSI-DATA直接对原始响应构造 Hankel 矩阵并做投影分解绕开了先算协方差矩阵的中间步骤抗噪能力和数值稳定性都更好。这篇围绕标题里的三个关键词展开数据驱动、随机子空间以及真正决定工程落地效果的模态自动识别——从代码到参数把整套流程说透。2. 数据驱动随机子空间的理论骨架Hankel矩阵、投影与SVD截断2.1 为什么选数据驱动SSI而不是协方差驱动SSI-COV做 SSI 之前先要选边站数据驱动还是协方差驱动。SSI-COV 先估计输出协方差矩阵再组装 Toeplitz 矩阵并做 SVDSSI-DATA 跳过协方差直接对 Hankel 矩阵做 QR 分解和投影。协方差这一步看着无害但协方差估计要把两两通道的乘积做时间平均信噪比低时协方差矩阵的秩会被测量噪声污染截断阶数难判断。数据驱动版本在 QR 分解阶段就让未来输出中与过去输出无关的噪声部分正交掉了投影矩阵的有效秩干净得多。现场数据里出现低频漂移、通道噪声不均匀时SSI-DATA 的稳定性优势尤其明显所以工程上我默认走数据驱动路线。SSI 的物理假设是结构响应可以写成一个离散状态空间模型x_{k1} A x_k w_ky_k C x_k v_k。w 是环境激励风、地脉动、交通荷载v 是测量噪声两者都被假设为零均值白噪声。随机子空间里的“随机”二字指的就是激励无法测量、只能用统计特性处理。A 矩阵的特征值决定系统固有频率和阻尼比C 矩阵与 A 的特征向量共同决定振型。整个识别过程的落点就是从输出矩阵 Y [y_0, y_1, …, y_{N-1}] 出发把 A 和 C 估出来。2.2 Hankel矩阵与QR分解把“过去”投影到“未来”Hankel 矩阵是 SSI 的起手式。核心思想是把输出信号按时间错位排列让行代表时间平移后的观测。常见做法是取 2i 个行块每个行块放全部 l 个通道前 i 个行块叫“过去”Y_p后 i 个行块叫“未来”Y_f整个矩阵尺寸为 2il × (N-2i1)。对这个大矩阵做 QR 分解后R 矩阵里 R21 那一块记录了“未来”中能被“过去”解释的部分这就是投影的本质——把未来的行空间投影到过去的行空间上。结构模态既然是系统内在属性它既影响过去也影响未来噪声没有这种跨时间的相关性会在投影里被压小。这个对比是 SSI 能在低信噪比下工作的根本原因。为什么用 QR 而不是最小二乘直接解QR 分解在数值上最稳定几百行乘几千列的矩阵在 numpy 或 MATLAB 里都是成熟优化过的内存占用和速度都可控。实际实现时对 [Y_p; Y_f]^T 做经济型 QR取 R 的右上块 R21 即可。2.3 SVD截断与系统阶数的第一道线索投影矩阵 O_i 的 SVD 分解 O_i U S V^T 是整个算法里信息量最大的一步。奇异值 σ_k 的衰减速度直接反映系统的有效阶数真实结构模态对应较大的奇异值噪声被压到后面的小奇异值上。import numpy as np # 假设投影矩阵 O 已经由 QR 分解构造完成 U, S, Vh np.linalg.svd(O, full_matricesFalse) # 看奇异值从第几个开始变平结构模态贡献大、拉高前几项 # 如果 S 在第 8 项后趋于平台说明系统本质上是 4 阶两自由度 print(S / S[0]) # 归一化奇异值谱肉眼找“悬崖”这段代码不截断任何东西只是把奇异值谱打出来。归一化之后结构模态对应的项通常在 1e-1 量级以上噪声平台在 1e-3 以下。两者之间如果有明显台阶阶数就好定了。但实桥上这个谱往往没有清晰悬崖尤其激励频带窄或传感器布点少时。所以阶数不能只靠 SVD 定稳定图才是最终裁决者——对应到标题里的“模态自动识别”这部分第四章会展开。3. 用Python实现SSI-DATA从仿真信号到频率阻尼振型3.1 生成一个有标准答案的仿真数据先造一个两自由度系统让识别结果有对照答案。用质量-弹簧-阻尼模型参数设定后理论频率约 3.1 Hz 和 8.1 Hz阻尼比在 1% 到 3% 之间。给质量 1 加白噪声激励两个质量的位移就是输出信号再叠一点测量噪声模拟现场。import numpy as np fs 50.0 # 采样率 50 Hz看 10 Hz 以内的模态足够 N 6000 # 120 秒数据 dt 1.0 / fs M np.eye(2) K np.array([[2000.0, -1000.0], [-1000.0, 1000.0]]) C 0.3 * M 0.0003 * K # Rayleigh 阻尼阻尼比约 1%~3% # 连续状态空间矩阵 A_c状态取 [位移, 速度] A_c np.vstack([ np.hstack([np.zeros((2, 2)), np.eye(2)]), np.hstack([-np.linalg.solve(M, K), -np.linalg.solve(M, C)]) ]) B np.array([0.0, 0.0, 1.0, 0.0]) # 激励作用在质量 1 上 # 白噪声激励 零阶保持离散化 A_d np.eye(4) A_c * dt state np.zeros(4) y np.zeros((2, N)) for k in range(N): w np.random.randn() state A_d state B * w * dt y[:, k] state[:2] 0.01 * np.random.randn(2)这里用零阶保持离散化近似dt 很小所以精度足够。叠加的 0.01 测量噪声对应约 1% 噪声水平比现场干净但足够用来验证算法主流程。输出 y 是 2×6000 的矩阵两行代表两个测点。3.2 SSI-DATA主流程实现核心函数按四步走构造 Hankel 矩阵、对合并矩阵做 QR、从 R21 得到投影并 SVD、截断后估 A 和 C。代码保持最小依赖只用 numpy 和 scipy。from scipy import linalg def ssi_data(y, i, n): 数据驱动随机子空间识别 y : 2D array, (测点数 l, 采样点数 T) i : Hankel 行块数每个行块含全部 l 个通道 n : SVD 截断阶数状态维数 返回状态矩阵 A_hat、输出矩阵 C_hat、奇异值谱 S l, T y.shape cols T - 2 * i 1 # 1. 构造 Hankel 矩阵2i 个行块第 k 块起始列是 k blocks [y[:, k:k cols] for k in range(2 * i)] Y_p np.vstack(blocks[:i]) # 前 i 块过去 Y_f np.vstack(blocks[i:]) # 后 i 块未来 # 2. 合并后做经济型 QR 分解取 R 的右上块 R21 Q, R linalg.qr(np.vstack([Y_p, Y_f]).T, modeeconomic) R21 R[i * l:, :i * l] # 3. 投影矩阵 O R21 Q1^T再 SVD O R21 Q[:, :i * l].T U, S, Vh linalg.svd(O, full_matricesFalse) # 4. 截断到 n 维构造扩展可观矩阵 Gamma U1 U[:, :n] S1 S[:n] Gamma U1 * np.sqrt(S1) # Gamma U1 * sqrt(S1) C_hat Gamma[:l, :] # 输出矩阵取 Gamma 前 l 行 X np.linalg.pinv(Gamma) Y_f # 卡尔曼滤波状态序列 # 用状态序列的时间错位最小二乘估计 A A_hat X[:, 1:] np.linalg.pinv(X[:, :-1]) return A_hat, C_hat, S关键在第 3 步对合并矩阵求 QRR21 是 Y_f 在 Y_p 行空间上的投影坐标再乘上 Q1 的转置才得到真正的投影矩阵。很多简化实现直接拿 R21 做 SVD严格说差了 Q1^T 这一项结果会有偏差。SVD 后截断到 n 维Gamma 的列数就是 n行数保持 il。C_hat 直接取 Gamma 的前 l 行因为扩展可观矩阵的第一块就是输出矩阵。3.3 状态矩阵A到模态参数log与特征向量估出 A_hat 和 C_hat 后频率、阻尼比、振型的提取是固定的三步操作def extract_modes(A_hat, C_hat, fs): # 1. 对状态矩阵做特征值分解 lam, psi np.linalg.eig(A_hat) # 2. 离散特征值转连续域s log(lam) * fs mu np.log(lam) * fs # 3. 频率和阻尼比 freq np.abs(mu) / (2 * np.pi) damp -mu.real / np.abs(mu) # 4. 振型C 矩阵映射到观测自由度 modes C_hat psi return freq, damp, modes离散状态矩阵的特征值 λ 和连续域特征值 s 之间的关系是 λ exp(s / fs)所以反解时 s log(λ)·fs。这里的 fs 单位要是 Hz别用角频率否则频率结果差 2π。阻尼比是连续特征值实部绝对值除以模长符号上注意负实部对应衰减。用两自由度仿真数据跑一遍设定 i20、n4得到的频率约 3.08 Hz 和 8.13 Hz阻尼比约 2% 和 1.5%与理论值吻合。第一次跑通后可以把 y 换成实测数据但这时候问题来了实测结构的真实阶数 n 是未知的不能拍脑袋设 4。这就是下一章“模态自动识别”要解决的。4. 模态自动识别稳定图、容差准则与层次聚类4.1 稳定图上的假模态从哪来n 取不同值时同样的数据会算出不同数量的“模态”。真实物理模态应该与 n 无关——不管系统阶数取 6 还是取 20第一阶频率都该出现在 3.1 Hz 附近。伪模态则随 n 变化跑来跑去位置不固定。把每个 n 对应的频率值画在横轴为频率、纵轴为 n 的图上真实模态形成一串垂直的“柱子”这就是稳定图的原理。稳定图上之所以总是混着散点主要三个原因激励并非理想白噪声风有低频着色、传感器噪声带入了非结构成分、以及非线性响应产生的谐波。自动识别要做的事就是把这些捣乱的散点和稳定的柱子分离开。4.2 稳定点判定频率、阻尼、MAC三参数容差表稳定图中“稳定”不是人眼说了算需要一组量化准则。参数容差物理含义频率变化 Δf 1%相对值固有频率是结构固有属性随阶数变化应该极小阻尼比变化 Δζ 5%相对值阻尼比本身难辨识容差放宽避免把真实模态滤掉振型 MAC 0.95两阶振型的空间相关性越低越可能是不同模态逐条解释一下。频率容差最严因为频率辨识最准1% 是常用起点模态密集时收紧到 0.5%。阻尼比识别方差大5% 相对容差对应工程经验太松会放进伪模态太紧会漏掉真实模态。MACModal Assurance Criterion是振型之间的相关系数0.95 以上视为同一个模态在相邻阶次间的稳定重现。三个条件同时满足才把一个点判为稳定点。4.3 用凝聚层次聚类自动挑出物理模态有了稳定点自动识别就变成一个聚类问题。把各阶次判定出的模态按频率作为特征做凝聚层次聚类频率接近的稳定点聚成一簇簇内点够多、散布够小就认为是物理模态。这里给一个可直接运行的版本。import numpy as np from scipy.cluster.hierarchy import linkage, fcluster def auto_select(freqs, damps, min_points5, tol_hz0.5): freqs / damps : 所有阶次下候选模态的频率与阻尼 min_points : 簇内最少点数少于该值视为伪模态 tol_hz : 聚类距离阈值单位 Hz X freqs.reshape(-1, 1) Z linkage(X, methodward) labels fcluster(Z, ttol_hz, criteriondistance) results [] for lab in sorted(set(labels)): idx labels lab if idx.sum() min_points: continue # 用簇内平均值作为该模态的最终估计 f_mean np.mean(freqs[idx]) d_mean np.mean(damps[idx]) results.append((f_mean, d_mean)) return results聚类特征用频率一维就够原因是稳定图柱子最显著的特征就是“同一频率附近重复出现”。ward 距离在这种一维小样本上表现稳定tol_hz 取 0.5 是经验值模态密集时调到 0.3结构稀疏时放宽到 1.0。min_points 设 5 意味着至少要在 5 个阶次上稳定出现低于这个数大多是一次性伪模态。阻尼比在聚类后取平均。这里有一个常见误区拿阻尼比当聚类特征。阻尼比估计离散度大同一个真实模态的阻尼可能从 1% 到 4% 分布当作聚类输入反而把簇打散。工程上我一般只用频率聚类再用阻尼比做事后过滤——如果某个“物理模态”的阻尼比超过 8%优先怀疑它是噪声着色形成的伪模态。5. 参数边界与现场排错让算法不给你“好看但错误”的结果5.1 采样率、数据长度与Hankel行数的配合参数之间的配合关系比单独调任何一个参数都重要。Hankel 矩阵的列数是 N - 2i 1行数是 2il。要保证 QR 分解和投影有意义列数至少要远大于行数经验上 N 4il 才算安全。用公式算6000 点、l2、i20 时行数是 80列数是 5961安全。如果 i 加到 80列数变成 5841行数 320仍然安全但再往上就危险了。采样率选择看目标频带。SSI 没有抗混叠滤波传感器采集时最好已经有了硬件滤波否则要先用数字滤波把高于关心的频率切掉。经验法则fs 至少是最高关注频率的 5 倍。50 Hz 采样看 10 Hz 以内的模态够用要看 30 Hz 就得提到 200 Hz 以上。数据长度则决定了统计置信度N 至少要覆盖最低关注模态的 10 个完整周期比如关注 0.5 Hz 模态就至少需要 20 秒数据实际建议留三倍余量。5.2 三个最常见的伪模态来源第一个来源是激励非白。环境激励的频谱不是纯平的比如风荷载在低频段能量高这会被 SSI 当成一个“模态”识别出来。特征是频率很低、阻尼比却很高、任意两个测点的振型分量同相——物理结构低频模态通常伴随明显振型变形不会所有测点一个方向。应对办法是提高阻尼比过滤阈值或者对比不同时段的识别结果物理模态位置不变激励谱产生的伪模态会漂移。第二个来源是谐波干扰。旋转机械的转频、电网 50 Hz 工频、桥梁涡振的锁定频率都会在稳定图上形成一条“永不消失的柱子”。谐波的典型特征是阻尼比极低往往低于 0.1%因为正弦信号没有衰减特性。用 damp 0.1% 做下限过滤能滤掉绝大部分谐波伪模态。第三个来源是传感器局部共振。加速度计安装不牢固时传感器本身在某个高频处振荡这个频率也是稳定的但它在不同测点之间没有一致的振型关系。用 MAC 验证如果某一阶模态在两个相邻测点之间的振型比值异常大先检查安装别急着改算法参数。5.3 自动识别结果可疑时的排查顺序结果不对时我一般按固定顺序排查。第一步看奇异值谱如果 S/S[0] 从第一项开始就滑得很平、没有明显台阶说明激励能量不够或测点信噪比太差这时候调聚类阈值没意义先解决数据质量。第二步看稳定图对比设置不同 tol_hz 时模态簇的数量变化数量稳定才可信。第三步核对 MAC把聚类得到的模态两两算 MAC如果不同频率的两个模态 MAC 超过 0.8说明振型信息不足测点数太少或测点位置靠近振型节点。网上流传的 SSI 速通资料不少链接已经打不开了与其在浏览器里反复刷不如把本文的仿真数据先跑通再替换成实测数据。自己亲手调一遍 i、n、tol_hz比收藏十个链接都管用这是五个常用参数里最值得花时间的一项。6. 验证技巧用MAC和交叉验证给识别结果上保险6.1 MAC振型一致性的量化依据MAC 是模态验证的通用语言公式为 MAC(a,b) |a^H b|² / (|a|²·|b|²)物理含义是两组振型向量的空间相关系数。同一模态在不同阶次下的 MAC 应接近 1不同模态间应明显低于 0.3。可以把 MAC 计算直接嵌入第四章的自动识别流程给聚类结果加一道事后闸门def mac(phi1, phi2): num np.abs(np.vdot(phi1, phi2)) den np.linalg.norm(phi1) * np.linalg.norm(phi2) return num / den # 聚类得到模态后把频率最接近两阶的振型做比较 # MAC 0.9 说明这两个识别结果确实是同一个模态 # MAC 0.7 说明测点信息不足需要复核布点MAC 阈值的使用场景要分清楚识别阶段用 0.95 滤稳定点验证阶段用 0.9 确认模态身份。后者更宽松因为验证面对的是已经筛选过的结果。6.2 留一段数据做交叉验证最实用的验证方法是把数据切两半分别独立跑一遍完整流程比较两段结果的频率差异。1 小时数据就前 30 分钟和后 30 分钟各跑一次频率偏差小于 0.5% 的模态可以信任偏差大于 2% 的基本判定为伪模态或非平稳影响。阻尼比的交叉验证容差放宽到 20%因为阻尼辨识方差本来就大。这种方法不增加任何硬件成本却能把伪模态漏网的概率降低一个数量级适用于从桥梁到风机叶片的各种尺度。6.3 降阶稳定图的技巧比纯频率聚类更稳的进阶做法是把稳定点画在“频率-阻尼”二维平面上然后对二维坐标做聚类。物理模态在二维图上聚成紧凑的点团伪模态则散布在平面各处比一维频率聚类更容易区分。代价是需要手动搜一遍二维聚类参数。实际操作里可以拿频率聚类的结果做初筛再用阻尼均值做核对如果某个簇内阻尼的标准差超过均值的一半果断删除——同一个物理模态在不同阶次上的阻尼估计可以波动但不该发散到这种程度。这套流程跑通之后SSI 的模态自动识别就能从“跑出图来人工挑”升级成“集群自动出报告”也就能长期稳定地铺到在线监测里去了。本文还有配套的精品资源点击获取 数据驱动随机子空间识别SSI：从Hankel矩阵到模态自动识别
简介面向结构动力学与健康监测研究者的数据驱动随机子空间识别DDRSI源码包聚焦模态自动识别可从含噪声的实测响应中直接提取频率、阻尼比和振型无需建立先验理论模型适用于桥梁、建筑、航空航天器等基础设施的振动监测与参数辨识。压缩包共24个文件包括13个MAT数据文件与11个M脚本数据文件涵盖多组仿真与实测信号如ASCE基准模型、连续6自由度结构等脚本则覆盖数据分块矩阵构建、稳定图绘制、模态参数自动筛选等关键环节便于对照学习和二次开发。已有1116人学习下载。资源包含完整可运行的主程序与配套函数并附有详细算例可帮助使用者系统掌握数据驱动随机子空间的建模步骤、去噪策略与虚假模态剔除方法为实际工程结构的状态评估和健康监测提供可复用的代码基础。1. 随机子空间识别不激励、不赋初值从实测响应里把模态拆出来一台风机机舱上的加速度计连续记录了几天环境振动数据没人敲它、也没开激振器却要报出结构前三阶频率、阻尼比和振型——这是运行模态分析里的典型任务。随机子空间识别SSI是这类任务里用得最广的算法族之一而数据驱动版本SSI-DATA直接对原始响应构造 Hankel 矩阵并做投影分解绕开了先算协方差矩阵的中间步骤抗噪能力和数值稳定性都更好。这篇围绕标题里的三个关键词展开数据驱动、随机子空间以及真正决定工程落地效果的模态自动识别——从代码到参数把整套流程说透。2. 数据驱动随机子空间的理论骨架Hankel矩阵、投影与SVD截断2.1 为什么选数据驱动SSI而不是协方差驱动SSI-COV做 SSI 之前先要选边站数据驱动还是协方差驱动。SSI-COV 先估计输出协方差矩阵再组装 Toeplitz 矩阵并做 SVDSSI-DATA 跳过协方差直接对 Hankel 矩阵做 QR 分解和投影。协方差这一步看着无害但协方差估计要把两两通道的乘积做时间平均信噪比低时协方差矩阵的秩会被测量噪声污染截断阶数难判断。数据驱动版本在 QR 分解阶段就让未来输出中与过去输出无关的噪声部分正交掉了投影矩阵的有效秩干净得多。现场数据里出现低频漂移、通道噪声不均匀时SSI-DATA 的稳定性优势尤其明显所以工程上我默认走数据驱动路线。SSI 的物理假设是结构响应可以写成一个离散状态空间模型x_{k1} A x_k w_ky_k C x_k v_k。w 是环境激励风、地脉动、交通荷载v 是测量噪声两者都被假设为零均值白噪声。随机子空间里的“随机”二字指的就是激励无法测量、只能用统计特性处理。A 矩阵的特征值决定系统固有频率和阻尼比C 矩阵与 A 的特征向量共同决定振型。整个识别过程的落点就是从输出矩阵 Y [y_0, y_1, …, y_{N-1}] 出发把 A 和 C 估出来。2.2 Hankel矩阵与QR分解把“过去”投影到“未来”Hankel 矩阵是 SSI 的起手式。核心思想是把输出信号按时间错位排列让行代表时间平移后的观测。常见做法是取 2i 个行块每个行块放全部 l 个通道前 i 个行块叫“过去”Y_p后 i 个行块叫“未来”Y_f整个矩阵尺寸为 2il × (N-2i1)。对这个大矩阵做 QR 分解后R 矩阵里 R21 那一块记录了“未来”中能被“过去”解释的部分这就是投影的本质——把未来的行空间投影到过去的行空间上。结构模态既然是系统内在属性它既影响过去也影响未来噪声没有这种跨时间的相关性会在投影里被压小。这个对比是 SSI 能在低信噪比下工作的根本原因。为什么用 QR 而不是最小二乘直接解QR 分解在数值上最稳定几百行乘几千列的矩阵在 numpy 或 MATLAB 里都是成熟优化过的内存占用和速度都可控。实际实现时对 [Y_p; Y_f]^T 做经济型 QR取 R 的右上块 R21 即可。2.3 SVD截断与系统阶数的第一道线索投影矩阵 O_i 的 SVD 分解 O_i U S V^T 是整个算法里信息量最大的一步。奇异值 σ_k 的衰减速度直接反映系统的有效阶数真实结构模态对应较大的奇异值噪声被压到后面的小奇异值上。import numpy as np # 假设投影矩阵 O 已经由 QR 分解构造完成 U, S, Vh np.linalg.svd(O, full_matricesFalse) # 看奇异值从第几个开始变平结构模态贡献大、拉高前几项 # 如果 S 在第 8 项后趋于平台说明系统本质上是 4 阶两自由度 print(S / S[0]) # 归一化奇异值谱肉眼找“悬崖”这段代码不截断任何东西只是把奇异值谱打出来。归一化之后结构模态对应的项通常在 1e-1 量级以上噪声平台在 1e-3 以下。两者之间如果有明显台阶阶数就好定了。但实桥上这个谱往往没有清晰悬崖尤其激励频带窄或传感器布点少时。所以阶数不能只靠 SVD 定稳定图才是最终裁决者——对应到标题里的“模态自动识别”这部分第四章会展开。3. 用Python实现SSI-DATA从仿真信号到频率阻尼振型3.1 生成一个有标准答案的仿真数据先造一个两自由度系统让识别结果有对照答案。用质量-弹簧-阻尼模型参数设定后理论频率约 3.1 Hz 和 8.1 Hz阻尼比在 1% 到 3% 之间。给质量 1 加白噪声激励两个质量的位移就是输出信号再叠一点测量噪声模拟现场。import numpy as np fs 50.0 # 采样率 50 Hz看 10 Hz 以内的模态足够 N 6000 # 120 秒数据 dt 1.0 / fs M np.eye(2) K np.array([[2000.0, -1000.0], [-1000.0, 1000.0]]) C 0.3 * M 0.0003 * K # Rayleigh 阻尼阻尼比约 1%~3% # 连续状态空间矩阵 A_c状态取 [位移, 速度] A_c np.vstack([ np.hstack([np.zeros((2, 2)), np.eye(2)]), np.hstack([-np.linalg.solve(M, K), -np.linalg.solve(M, C)]) ]) B np.array([0.0, 0.0, 1.0, 0.0]) # 激励作用在质量 1 上 # 白噪声激励 零阶保持离散化 A_d np.eye(4) A_c * dt state np.zeros(4) y np.zeros((2, N)) for k in range(N): w np.random.randn() state A_d state B * w * dt y[:, k] state[:2] 0.01 * np.random.randn(2)这里用零阶保持离散化近似dt 很小所以精度足够。叠加的 0.01 测量噪声对应约 1% 噪声水平比现场干净但足够用来验证算法主流程。输出 y 是 2×6000 的矩阵两行代表两个测点。3.2 SSI-DATA主流程实现核心函数按四步走构造 Hankel 矩阵、对合并矩阵做 QR、从 R21 得到投影并 SVD、截断后估 A 和 C。代码保持最小依赖只用 numpy 和 scipy。from scipy import linalg def ssi_data(y, i, n): 数据驱动随机子空间识别 y : 2D array, (测点数 l, 采样点数 T) i : Hankel 行块数每个行块含全部 l 个通道 n : SVD 截断阶数状态维数 返回状态矩阵 A_hat、输出矩阵 C_hat、奇异值谱 S l, T y.shape cols T - 2 * i 1 # 1. 构造 Hankel 矩阵2i 个行块第 k 块起始列是 k blocks [y[:, k:k cols] for k in range(2 * i)] Y_p np.vstack(blocks[:i]) # 前 i 块过去 Y_f np.vstack(blocks[i:]) # 后 i 块未来 # 2. 合并后做经济型 QR 分解取 R 的右上块 R21 Q, R linalg.qr(np.vstack([Y_p, Y_f]).T, modeeconomic) R21 R[i * l:, :i * l] # 3. 投影矩阵 O R21 Q1^T再 SVD O R21 Q[:, :i * l].T U, S, Vh linalg.svd(O, full_matricesFalse) # 4. 截断到 n 维构造扩展可观矩阵 Gamma U1 U[:, :n] S1 S[:n] Gamma U1 * np.sqrt(S1) # Gamma U1 * sqrt(S1) C_hat Gamma[:l, :] # 输出矩阵取 Gamma 前 l 行 X np.linalg.pinv(Gamma) Y_f # 卡尔曼滤波状态序列 # 用状态序列的时间错位最小二乘估计 A A_hat X[:, 1:] np.linalg.pinv(X[:, :-1]) return A_hat, C_hat, S关键在第 3 步对合并矩阵求 QRR21 是 Y_f 在 Y_p 行空间上的投影坐标再乘上 Q1 的转置才得到真正的投影矩阵。很多简化实现直接拿 R21 做 SVD严格说差了 Q1^T 这一项结果会有偏差。SVD 后截断到 n 维Gamma 的列数就是 n行数保持 il。C_hat 直接取 Gamma 的前 l 行因为扩展可观矩阵的第一块就是输出矩阵。3.3 状态矩阵A到模态参数log与特征向量估出 A_hat 和 C_hat 后频率、阻尼比、振型的提取是固定的三步操作def extract_modes(A_hat, C_hat, fs): # 1. 对状态矩阵做特征值分解 lam, psi np.linalg.eig(A_hat) # 2. 离散特征值转连续域s log(lam) * fs mu np.log(lam) * fs # 3. 频率和阻尼比 freq np.abs(mu) / (2 * np.pi) damp -mu.real / np.abs(mu) # 4. 振型C 矩阵映射到观测自由度 modes C_hat psi return freq, damp, modes离散状态矩阵的特征值 λ 和连续域特征值 s 之间的关系是 λ exp(s / fs)所以反解时 s log(λ)·fs。这里的 fs 单位要是 Hz别用角频率否则频率结果差 2π。阻尼比是连续特征值实部绝对值除以模长符号上注意负实部对应衰减。用两自由度仿真数据跑一遍设定 i20、n4得到的频率约 3.08 Hz 和 8.13 Hz阻尼比约 2% 和 1.5%与理论值吻合。第一次跑通后可以把 y 换成实测数据但这时候问题来了实测结构的真实阶数 n 是未知的不能拍脑袋设 4。这就是下一章“模态自动识别”要解决的。4. 模态自动识别稳定图、容差准则与层次聚类4.1 稳定图上的假模态从哪来n 取不同值时同样的数据会算出不同数量的“模态”。真实物理模态应该与 n 无关——不管系统阶数取 6 还是取 20第一阶频率都该出现在 3.1 Hz 附近。伪模态则随 n 变化跑来跑去位置不固定。把每个 n 对应的频率值画在横轴为频率、纵轴为 n 的图上真实模态形成一串垂直的“柱子”这就是稳定图的原理。稳定图上之所以总是混着散点主要三个原因激励并非理想白噪声风有低频着色、传感器噪声带入了非结构成分、以及非线性响应产生的谐波。自动识别要做的事就是把这些捣乱的散点和稳定的柱子分离开。4.2 稳定点判定频率、阻尼、MAC三参数容差表稳定图中“稳定”不是人眼说了算需要一组量化准则。参数容差物理含义频率变化 Δf 1%相对值固有频率是结构固有属性随阶数变化应该极小阻尼比变化 Δζ 5%相对值阻尼比本身难辨识容差放宽避免把真实模态滤掉振型 MAC 0.95两阶振型的空间相关性越低越可能是不同模态逐条解释一下。频率容差最严因为频率辨识最准1% 是常用起点模态密集时收紧到 0.5%。阻尼比识别方差大5% 相对容差对应工程经验太松会放进伪模态太紧会漏掉真实模态。MACModal Assurance Criterion是振型之间的相关系数0.95 以上视为同一个模态在相邻阶次间的稳定重现。三个条件同时满足才把一个点判为稳定点。4.3 用凝聚层次聚类自动挑出物理模态有了稳定点自动识别就变成一个聚类问题。把各阶次判定出的模态按频率作为特征做凝聚层次聚类频率接近的稳定点聚成一簇簇内点够多、散布够小就认为是物理模态。这里给一个可直接运行的版本。import numpy as np from scipy.cluster.hierarchy import linkage, fcluster def auto_select(freqs, damps, min_points5, tol_hz0.5): freqs / damps : 所有阶次下候选模态的频率与阻尼 min_points : 簇内最少点数少于该值视为伪模态 tol_hz : 聚类距离阈值单位 Hz X freqs.reshape(-1, 1) Z linkage(X, methodward) labels fcluster(Z, ttol_hz, criteriondistance) results [] for lab in sorted(set(labels)): idx labels lab if idx.sum() min_points: continue # 用簇内平均值作为该模态的最终估计 f_mean np.mean(freqs[idx]) d_mean np.mean(damps[idx]) results.append((f_mean, d_mean)) return results聚类特征用频率一维就够原因是稳定图柱子最显著的特征就是“同一频率附近重复出现”。ward 距离在这种一维小样本上表现稳定tol_hz 取 0.5 是经验值模态密集时调到 0.3结构稀疏时放宽到 1.0。min_points 设 5 意味着至少要在 5 个阶次上稳定出现低于这个数大多是一次性伪模态。阻尼比在聚类后取平均。这里有一个常见误区拿阻尼比当聚类特征。阻尼比估计离散度大同一个真实模态的阻尼可能从 1% 到 4% 分布当作聚类输入反而把簇打散。工程上我一般只用频率聚类再用阻尼比做事后过滤——如果某个“物理模态”的阻尼比超过 8%优先怀疑它是噪声着色形成的伪模态。5. 参数边界与现场排错让算法不给你“好看但错误”的结果5.1 采样率、数据长度与Hankel行数的配合参数之间的配合关系比单独调任何一个参数都重要。Hankel 矩阵的列数是 N - 2i 1行数是 2il。要保证 QR 分解和投影有意义列数至少要远大于行数经验上 N 4il 才算安全。用公式算6000 点、l2、i20 时行数是 80列数是 5961安全。如果 i 加到 80列数变成 5841行数 320仍然安全但再往上就危险了。采样率选择看目标频带。SSI 没有抗混叠滤波传感器采集时最好已经有了硬件滤波否则要先用数字滤波把高于关心的频率切掉。经验法则fs 至少是最高关注频率的 5 倍。50 Hz 采样看 10 Hz 以内的模态够用要看 30 Hz 就得提到 200 Hz 以上。数据长度则决定了统计置信度N 至少要覆盖最低关注模态的 10 个完整周期比如关注 0.5 Hz 模态就至少需要 20 秒数据实际建议留三倍余量。5.2 三个最常见的伪模态来源第一个来源是激励非白。环境激励的频谱不是纯平的比如风荷载在低频段能量高这会被 SSI 当成一个“模态”识别出来。特征是频率很低、阻尼比却很高、任意两个测点的振型分量同相——物理结构低频模态通常伴随明显振型变形不会所有测点一个方向。应对办法是提高阻尼比过滤阈值或者对比不同时段的识别结果物理模态位置不变激励谱产生的伪模态会漂移。第二个来源是谐波干扰。旋转机械的转频、电网 50 Hz 工频、桥梁涡振的锁定频率都会在稳定图上形成一条“永不消失的柱子”。谐波的典型特征是阻尼比极低往往低于 0.1%因为正弦信号没有衰减特性。用 damp 0.1% 做下限过滤能滤掉绝大部分谐波伪模态。第三个来源是传感器局部共振。加速度计安装不牢固时传感器本身在某个高频处振荡这个频率也是稳定的但它在不同测点之间没有一致的振型关系。用 MAC 验证如果某一阶模态在两个相邻测点之间的振型比值异常大先检查安装别急着改算法参数。5.3 自动识别结果可疑时的排查顺序结果不对时我一般按固定顺序排查。第一步看奇异值谱如果 S/S[0] 从第一项开始就滑得很平、没有明显台阶说明激励能量不够或测点信噪比太差这时候调聚类阈值没意义先解决数据质量。第二步看稳定图对比设置不同 tol_hz 时模态簇的数量变化数量稳定才可信。第三步核对 MAC把聚类得到的模态两两算 MAC如果不同频率的两个模态 MAC 超过 0.8说明振型信息不足测点数太少或测点位置靠近振型节点。网上流传的 SSI 速通资料不少链接已经打不开了与其在浏览器里反复刷不如把本文的仿真数据先跑通再替换成实测数据。自己亲手调一遍 i、n、tol_hz比收藏十个链接都管用这是五个常用参数里最值得花时间的一项。6. 验证技巧用MAC和交叉验证给识别结果上保险6.1 MAC振型一致性的量化依据MAC 是模态验证的通用语言公式为 MAC(a,b) |a^H b|² / (|a|²·|b|²)物理含义是两组振型向量的空间相关系数。同一模态在不同阶次下的 MAC 应接近 1不同模态间应明显低于 0.3。可以把 MAC 计算直接嵌入第四章的自动识别流程给聚类结果加一道事后闸门def mac(phi1, phi2): num np.abs(np.vdot(phi1, phi2)) den np.linalg.norm(phi1) * np.linalg.norm(phi2) return num / den # 聚类得到模态后把频率最接近两阶的振型做比较 # MAC 0.9 说明这两个识别结果确实是同一个模态 # MAC 0.7 说明测点信息不足需要复核布点MAC 阈值的使用场景要分清楚识别阶段用 0.95 滤稳定点验证阶段用 0.9 确认模态身份。后者更宽松因为验证面对的是已经筛选过的结果。6.2 留一段数据做交叉验证最实用的验证方法是把数据切两半分别独立跑一遍完整流程比较两段结果的频率差异。1 小时数据就前 30 分钟和后 30 分钟各跑一次频率偏差小于 0.5% 的模态可以信任偏差大于 2% 的基本判定为伪模态或非平稳影响。阻尼比的交叉验证容差放宽到 20%因为阻尼辨识方差本来就大。这种方法不增加任何硬件成本却能把伪模态漏网的概率降低一个数量级适用于从桥梁到风机叶片的各种尺度。6.3 降阶稳定图的技巧比纯频率聚类更稳的进阶做法是把稳定点画在“频率-阻尼”二维平面上然后对二维坐标做聚类。物理模态在二维图上聚成紧凑的点团伪模态则散布在平面各处比一维频率聚类更容易区分。代价是需要手动搜一遍二维聚类参数。实际操作里可以拿频率聚类的结果做初筛再用阻尼均值做核对如果某个簇内阻尼的标准差超过均值的一半果断删除——同一个物理模态在不同阶次上的阻尼估计可以波动但不该发散到这种程度。这套流程跑通之后SSI 的模态自动识别就能从“跑出图来人工挑”升级成“集群自动出报告”也就能长期稳定地铺到在线监测里去了。本文还有配套的精品资源点击获取