简介:这是一份基于随机子空间与稳定性选择的MATLAB工具包,面向高维数据特征筛选、模型泛化性分析与机器学习建模场景,通过多次随机抽取特征子空间训练模型,计算特征在不同子集中的出现频率,从而绘制出直观的SSI稳定图,帮助研究者识别对结果贡献最稳定的特征。压缩包共四个m脚本,整体大小约四KB,涵盖主运行程序、数据滤波与降采样、协方差计算等模块,结构紧凑;已有四百零八人学习或下载,适合数据预处理、特征工程和信号分析方向的开发人员参考。解压后运行主脚本即可启动完整流程:从加载数据、随机子空间生成,到模型训练、稳定性打分与绘图,一步完成;同时附带的预处理脚本可对原始信号进行降采样和滤波,减少噪声干扰,提升后续建模的稳健性与解释力。该工具在维度高、样本规模适中的场景下尤为实用,能显著降低全空间遍历带来的计算开销,并为后续建模提供明确的有效特征清单。
1. SSI稳定图到底在解决什么问题:从一堆振动波形里逼出真实模态
拿到一份标着SSI.rar的实测数据包时,我第一个动作不是解压找代码,而是先问一句:这里面的随机子空间识别脚本,到底能不能输出一张稳定图?因为对做实测模态分析的人来说,SSI(随机子空间识别)和稳定图是绑在一起的。SSI负责从环境激励下的加速度响应里建状态空间模型,稳定图负责回答“哪些识别结果是可信的模态”。没有稳定图的SSI,就像带回家一堆滤波器参数却不知道哪个频点代表结构自身频率,只能靠猜。现实场景往往是:桥梁、风机塔筒或高层建筑长期监测,只有输出信号没有输入力,频域峰值法抓出来的谱峰混杂着电磁干扰和倍频谐波;而稳定图上,真实模态会在不同系统阶次下反复出现,连成一条竖直的稳定柱。这篇文章就从原理、可复现代码、参数标定到避坑,把这套流程完整拆给你。
2. 随机子空间识别的数学骨架:状态空间、Toeplitz矩阵与奇异值截断
2.1 为什么偏偏选随机子空间而不是频域法
实测模态分析里,传统频域法比如峰值拾取、最小二乘复频域法,依赖频响函数或输出谱。可环境激励下没有可测的输入力,只能把响应谱当作频响函数的替代,这时峰值周围往往被噪声抹平,阻尼比估值也很不稳定。更麻烦的是,密集模态在谱线上粘连成一个小包,峰值法根本数不过去。SSI是纯时域方法,它把未知的环境激励统一建模为白噪声过程噪声,不需要输入力,直接从输出的二阶统计量里重构系统状态矩阵。随机子空间这个名字听着玄,实际做的事就是把响应信号投影到由过去输出张成的子空间里,用奇异值分解把一个高维随机过程压缩成低阶确定性的状态空间模型。它不像频域法那样需要人为选频带,用一段足够长的加速度时程就能同时给出频率、阻尼比和振型。对桥梁、风电塔、大跨结构这类环境激励占主导的场景,SSI几乎是默认选项。
当然它也有代价:阶次难定、对数据长度敏感、稳定图需要判读。所以业界才把“系统辨识”和“稳定图”绑定在一起,用后者的重复性来弥补前者的不确定性。理解了这一点,你就知道为什么每个SSI工具包里,稳定图脚本往往比识别算法本体还长。
2.2 状态空间模型与协方差序列构建
随机子空间识别的起点是离散状态空间方程:
x_{k+1} = A x_k + w_k
y_k = C x_k + v_k
其中 x_k 是 n 维状态向量,y_k 是 m 维观测输出(比如多个加速度测点),w_k 是过程噪声,v_k 是测量噪声,两者都假设为零均值白噪声。环境激励的合力被吸收进 w_k,所以模型不需要显式输入项。我们的目标就是从输出序列 y_k 里估计状态矩阵 A 和输出矩阵 C。结构模态信息全部藏在 A 的特征值里:设离散特征值为 λ,采样间隔为 dt,则连续域极点 s = log(λ)/dt,模态频率 f = |s|/(2π),阻尼比 ζ = -Re(s)/|s|。这个映射在代码里只有两行,但它是整个识别的终点。
要从随机响应里还原 A,需要利用输出序列的协方差。定义时滞 l 的协方差矩阵 R_l = E[y_{k+l} y_k^T],实际工程中用有限样本估计。协方差序列 R_1, R_2, ... 携带了系统的脉冲响应信息,但直接看单个 R_l 很难提取模态。常见做法是把它们堆成一个块 Toeplitz 矩阵 T,块的排列方式近似于将不同时滞的两两相关打包在一起。对 T 做奇异值分解,就能找到“可观测子空间”。
块 Toeplitz 矩阵的具体形式:取块行数 i(一般取 10~50),构造维度为 (i·m)×(i·m) 的矩阵,第 p 行第 q 列的块取 R_{p+q}(这里 R 的下标按从 1 开始计)。也就是说 T 的每个块都是由 m×m 的滞后协方差矩阵填充。这样构造出来的矩阵把多个时滞的统计信息叠加,等价于把过去和未来输出做投影。和复杂的 Data-SSI 投影算法相比,这种协方差 SSI(Cov-SSI)形式实现简单、稳定性好,适合快速验证,也是很多开源工具包的底层。
2.3 SVD截断与阶次确定的实操逻辑
对块 Toeplitz 矩阵 T 做奇异值分解:T = U S V^T。理想情况下,系统真实阶次为 n_0 时,S 中只有 n_0 个非零奇异值。但实测数据总有噪声,奇异值尾部不会归零,而是缓慢衰减。要真正确定阶次,得靠稳定图扫描,而不是眼睛盯着奇异值的“悬崖”。不过,SVD 截断的思路仍然贯穿算法:取前 n 个奇异值对应的左右奇异向量,构造可观测性矩阵 Γ = U(:,1:n) · S(1:n)^(1/2)。然后从 Γ 中提取 C 和 A:C 直接取 Γ 的前 m 行,A 利用 Γ 的移不变性,用最小二乘解方程 Γ_up · A = Γ_down,其中 Γ_up 是 Γ 去掉最后 m 行,Γ_down 是去掉前 m 行。
这里有几个直接影响结果的操作参数:
| 参数 | 作用 | 经验范围 |
|---|---|---|
| 块行数 i | 决定 Toeplitz 矩阵的尺寸和可观测性范围 | 10~50,且 i 需大于最大目标阶次 |
| 扫描阶次 n | 每个 i 下从 1 递增到最大阶次,形成稳定图纵轴 | 1~2·i·m 之间,通常取 1~min(50, i·m/2) |
| 数据长度 N | 影响协方差统计精度 | 至少 1000~2000 点,建议响应循环数超过 100 个周期 |
| 采样率 fs | 确定奈奎斯特频率和离散模型带宽 | 覆盖目标最高模态频率的 5~10 倍以上 |
为什么用 SVD 而不是特征值分解?因为 Toeplitz 矩阵通常接近病态,特征值分解对噪声敏感;SVD 从左右奇异向量同时提取可观测子空间,数值稳定性好。在实际项目中,我一般先把最大扫描阶数设为 i·m/2,再观察稳定图里稳定点集中的阶数范围。如果最大阶数设得太低,真实高阶模态会被截掉;设得太高,系统矩阵 A 的维度变大,噪声制造的虚假极点会剧增,稳定图会变成“雪花屏”。
小节之间内容要连贯。这里想强调:随机子空间识别的核心不是“子空间”这个名词,而是用统计量构造可辨识的数据矩阵,再用稳定的线性代数工具提取状态矩阵。下一章,我们直接把上述理论转成可在本地运行的 Python 代码,并用一个已知两自由度系统做验证。
3. 用一段可复现的代码跑通SSI识别:从加速度时程到稳定图
3.1 构造一段可复现的模拟数据(两自由度系统)
写 SSI 代码之前,先准备一组“已知答案”的数据。模拟一个两自由度质量-弹簧-阻尼系统,真实频率设为 2.1 Hz 和 5.8 Hz,阻尼比分别 2% 和 3%。用随机白噪声激励两个自由度,采样率 100 Hz,时长 120 秒,生成两组加速度响应作为“实测输出”。这样跑完 SSI 后,我们能把识别结果和理论值对比,检验算法是否正确。
import numpy as np from scipy.signal import StateSpace, lsim # 两自由度系统参数 M = np.array([[2.0, 0.0], [0.0, 1.0]]) C = np.array([[3.0, -0.8], [-0.8, 1.5]]) K = np.array([[300.0, -120.0], [-120.0, 180.0]]) # 连续状态空间:状态 x = [位移1, 位移2, 速度1, 速度2] A_cont = np.block([[np.zeros((2, 2)), np.eye(2)], [-np.linalg.solve(M, K), -np.linalg.solve(M, C)]]) # 外力直接作用在两个自由度上,输入维度=2 B_cont = np.vstack([np.zeros((2, 2)), np.linalg.solve(M, np.eye(2))]) # 加速度输出,维度=2 C_cont = np.block([-np.linalg.solve(M, K), -np.linalg.solve(M, C)]) D_cont = np.linalg.solve(M, np.eye(2)) sys = StateSpace(A_cont, B_cont, C_cont, D_cont) dt = 0.01 t = np.arange(0, 120, dt) rng = np.random.default_rng(42) u = rng.normal(0, 1, (len(t), 2)) # 两个自由度的独立白噪声激励 _, _, y = lsim(sys, u, t) y = y.T # 转成 (2, N),便于按通道处理参数说明:A_cont是 4×4 连续状态矩阵,由质量、阻尼、刚度矩阵拼装而成;B_cont把外力映射到加速度方程;输出矩阵C_cont直接导出加速度,省去了从位移响应求导的麻烦。lsim用离散采样激励序列驱动连续系统,返回 (N, 2) 的响应,我们转置成 (2, N)。运行后,y[0]和y[1]就是两个测点的加速度时程。
为了确认数据可用,可以快速画频谱看峰值,但这里省略画图步骤。直接用下面的 SSI 函数跑识别,重点在验证算法是否还原出 2.1 Hz 和 5.8 Hz。
3.2 SSI核心识别函数:协方差Toeplitz与SVD提取
以下函数实现 Cov-SSI 的核心流程。输入加速度矩阵、块行数 i、当前阶次 n、采样间隔 dt,输出从该阶次模型里解出的模态频率、阻尼比和振型矩阵。
def cov_ssi(y, i, n, dt): """ y: (m, N) 输出时程,m为测点数 i: 块行数 n: 当前系统阶次 dt: 采样间隔,单位秒 返回: 频率数组(Hz),阻尼比数组(百分比),振型矩阵(m, n) """ m, N = y.shape # 1. 计算滞后协方差序列 R_1 ... R_{2i-1} max_lag = 2 * i - 1 R = [] for lag in range(1, max_lag + 1): y0 = y[:, :N-lag] - y[:, :N-lag].mean(axis=1, keepdims=True) y1 = y[:, lag:] - y[:, lag:].mean(axis=1, keepdims=True) R.append((y1 @ y0.T) / (N - lag)) # 2. 构建块Toeplitz矩阵 T (i*m, i*m) T = np.zeros((i * m, i * m)) for p in range(i): for q in range(i): idx = p + q if idx < len(R): T[p*m:(p+1)*m, q*m:(q+1)*m] = R[idx] # 3. SVD分解并截断到前n阶 U, S, _ = np.linalg.svd(T) U1 = U[:, :n] S1 = S[:n] Gamma = U1 @ np.diag(np.sqrt(S1)) # 4. 提取输出矩阵 C 和状态矩阵 A C_mat = Gamma[:m, :] Gamma_up = Gamma[:m*(i-1), :] Gamma_down = Gamma[m:, :] A_mat, _, _, _ = np.linalg.lstsq(Gamma_up, Gamma_down, rcond=None) # 5. 特征值分解,映射到连续域 eigvals, eigvecs = np.linalg.eig(A_mat) s = np.log(eigvals) / dt freq = np.abs(s) / (2 * np.pi) damp = -np.real(s) / np.abs(s) * 100 # 6. 振型 = C_mat * 特征向量,这里直接返回C_mat和特征向量的乘积 modes = C_mat @ eigvecs return freq, damp, modes逻辑说明:第 1 步计算滞后协方差时,先减去均值再算相关,避免直流分量污染。第 2 步的块索引idx = p + q对应 R_{p+q+1},例如 p=0, q=0 取 R_1,p=i-1, q=i-1 取 R_{2i-1},正好覆盖需要的滞时范围。第 3 步用 SVD 截断得到可观测性矩阵 Γ,这一步是随机子空间的核心,它把统计信息压缩成低秩状态空间。第 4 步用 Γ 的移位不变性求 A,最小二乘会被自动处理为超定方程。第 5 步注意np.log可能对接近实轴的复特征值产生噪动,但通常只取正频率部分。
参数选择上,i不能太小,否则 Toeplitz 矩阵提供的信息不足,i太大则计算量大且协方差估值噪声增加。对两个测点、最大阶次 40 的情况,取i=30足够。如果识别高阶模态,i要相应增大,并确保i*m大于n_max。
3.3 模态参数提取与稳定图数据生成
单次调用cov_ssi只能得到特定阶次下的模态“猜测”。稳定图的思路就是让阶次 n 从 1 逐步增大到最大值,每个阶次都做一次 SVD 和特征分解,再把结果收集起来。代码如下:
# 扫描阶次,生成稳定图原始数据 max_order = 40 i_blocks = 30 stability_data = [] for n in range(1, max_order + 1): freq, damp, modes = cov_ssi(y, i_blocks, n, dt) # 只保留正频率且落在0.1~10Hz范围内的真实候选 mask = (freq > 0.1) & (freq < 10.0) stability_data.append({ 'order': n, 'freq': freq[mask], 'damp': damp[mask], 'modes': modes[:, mask] }) # 把每个阶次的所有频率点打平,用于后续绘图/聚类 all_freq = np.concatenate([d['freq'] for d in stability_data])参数说明:max_order=40是指扫描最高系统阶次,理论占空比不超过i_blocks * 测点数 / 2 = 30*2/2=30,所以 40 已经超过可辨识上限,会引入较多噪声模态,但这正是稳定图需要的——噪声模态在阶次变化时会飘动,真实模态则稳定驻留。i_blocks=30是块行数,取的偏大是为了让 Toeplitz 矩阵能容纳 40 阶模型;块行数必须大于最大阶次,否则高阶截断失真。
稳定图原始数据并不是最终答案。如果把所有阶次的所有频率点都画出来,会看到密密麻麻的散点,完全无法区分真假。下一步必须用稳定判据把“相邻阶次间保持一致的模态”找出来,这才是稳定图的真正含义。下面进入第 4 章,把判据和容差一次说清。
4. 稳定图的参数标定与虚假模态剔除规则
4.1 三个稳定判据:频率、阻尼、MAC
稳定图上的每个点代表某一阶次下识别出的一个模态。如果两个连续阶次(通常是 n 和 n-1)识别出的模态频率几乎不变、阻尼比接近、振型高度相关,就称它为“稳定点”。工程判据有三个,缺一不可:
| 判据 | 常用阈值 | 说明 |
|---|---|---|
| 频率误差 | |f_n - f_{n-1}| / f_{n-1} ≤ 1% | 最核心、最稳定的判据,噪声大时可放宽到 2%~3% |
| 阻尼比误差 | |ζ_n - ζ_{n-1}| ≤ 5%(绝对值)或 ≤ 10%(相对值) | 阻尼比本身离散性大,只用它做辅助判据 |
| 振型相关 MAC | MAC ≥ 0.8~0.95 | 振型向量需要做归一化,传感器少时 MAC 容易被高估,需谨慎 |
为什么要同时用三个?因为频率稳定并不代表模态是真的。噪声主导的伪模态也可能在某个窄带内反复出现,但它的阻尼比会剧烈跳动,振型向量方向随机。只有频率、阻尼、振型同时稳定,才说明这个极点对应的状态空间模型在阶次扩展时被反复捕捉到,大概率是物理模态。
MAC 的计算公式为:MAC = |a^T b|^2 / ((a^T a)(b^T b)),其中 a、b 是同一阶次的振型向量。MAC 接近 1 表示两个振型几乎同向,接近 0 表示正交。实际代码里建议先对模态向量按某个测点归一化,再算 MAC,避免幅值影响。
4.2 容差怎么设:从经验值到自适应
容差设置没有绝对标准,我的经验是分两步:先用严格的默认容差生成稳定图,再根据稳定点的“柱高”(同一频率附近稳定点重复次数)调整。默认参数如下:
- 频率容差:0.01(即 1%)
- 阻尼容差:绝对值 5%
- MAC 容差:0.8
如果这样得到的稳定图几乎没有稳定点,通常是数据长度不足或噪声太大,这时把频率容差放宽到 0.02,阻尼容差放宽到 10%,MAC 降到 0.7。注意放宽不等于降低判据质量,而是换取识别覆盖率,最终以交叉验证结果为准。
实现稳定点判定的最小代码如下:
def mark_stable_points(stability_data, i_blocks, max_order, tol_freq=0.01, tol_damp=5.0, tol_mac=0.80): """ 输入stability_data(第3.3节的列表),返回每个阶次的标记矩阵 标记规则:当前阶次n与前一阶次n-1的模态匹配则稳定 """ m = 2 # 测点数,实际应从y.shape[0]读 stable_info = [] prev_freq = None prev_damp = None prev_modes = None for d in stability_data: cur_freq = d['freq'] cur_damp = d['damp'] cur_modes = d['modes'] marks = [] if prev_freq is None or len(cur_freq) == 0: marks = np.zeros(len(cur_freq), dtype=bool) else: for j, f in enumerate(cur_freq): # 在上一阶次找最近频率 err = np.abs(prev_freq - f) / f cand = np.where(err < tol_freq)[0] stable = False for c in cand: # 阻尼差 if abs(cur_damp[j] - prev_damp[c]) > tol_damp: continue # 振型MAC a = cur_modes[:, j].reshape(-1, 1) b = prev_modes[:, c].reshape(-1, 1) mac = np.abs(a.conj().T @ b)**2 / ((a.conj().T @ a) * (b.conj().T @ b)) if mac.item() >= tol_mac: stable = True break marks[j] = stable stable_info.append(marks) prev_freq, prev_damp, prev_modes = cur_freq, cur_damp, cur_modes return stable_info这段代码的核心是“与前一阶次匹配”。它有一个隐含假设:当阶次增加一阶时,真实模态的位置只会微动,而噪声模态可能跳到完全不同的频率。匹配时先按频率卡一道门槛,再用阻尼和 MAC 卡第二道门槛,三道全过才标稳定。实际应用中可以对比 n 和 n-2 甚至 n-3,但相邻阶次最简单,效果也够。
4.3 手动判读与自动聚类
稳定点标记完成后,画稳定图通常用横轴频率、纵轴阶次,稳定点用一个加粗的圆点标出。真实模态在图上表现为一条竖直的“柱子”:同一个频率附近,多个阶次都出现稳定点。手工判读时,我习惯先找柱子,再看柱子对应的阻尼是否落在合理范围(钢结构 0.3%~2%,混凝土结构 1%~5% 左右),最后从稳定柱所在阶次里取模态参数。
人工判读存在主观性,尤其面对几十个候选柱时。所以常用自动聚类做初筛。最简单有效的方法是频率直方图统计:把稳定点的频率按 0.02 Hz 的 bin 分箱,统计每个 bin 的稳定点个数,超过阈值(比如 3 个)就认为该频率附近可能存在真实模态。代码如下:
def auto_pick_frequencies(stable_points, all_freq_bin, min_hits=3, bin_width=0.02): """根据稳定点频率直方图粗选真实模态""" hits, bin_edges = np.histogram(all_freq_bin, bins=np.arange(0, 10, bin_width)) candidates = bin_edges[:-1][hits >= min_hits] # 合并相邻bin的中心频率(简单去重) merged = [] for f in candidates: if not merged or f - merged[-1] > bin_width * 1.5: merged.append(f) return merged参数min_hits=3表示至少 3 个阶次在该频率 bin 里产生稳定点,这个值越高越保守。自动聚类只能当“草稿”,最终还是要回到稳定图上人工核对,尤其是两个稳定柱距离很近时,要检查它们是不是同一阶的重复峰值。
稳定图的主观性不会完全消失,但把容差和聚类规则固定下来,同一批数据在不同人手里的结果就能基本对齐。这也是企业级模态分析报告里稳定图被反复强调的原因:它把黑匣子算法的输出变成了可审计的证据。
5. SSI实操避坑与常见问题排查
5.1 阻尼比系统性偏大的陷阱
现象:SSI 识别出的阻尼比明显大于数值模型或衰减法结果,比如钢桥识别出 3% 的阻尼,但理论只有 0.8%。
原因:环境激励不是理想白噪声,常含低频风压漂移和窄带附加振动;测量噪声会提高系统表观阻尼;此外带通滤波器如果设置不恰当,会在目标频带边缘引入相位畸变,让极点向复平面内压缩。
解决:在跑 SSI 之前对加速度时程做去趋势和零相位带通滤波,滤波带宽要比目标模态频率宽 1.5 倍以上;增加数据长度,让协方差估计更稳定;对比多段独立数据,若阻尼均值比单段结果低,说明单段阻尼被噪声抬高。还有一种常见做法是只保留稳定图中阻尼值相对稳定(相邻阶次阻尼差小)的模态,舍弃那些阻尼随阶次剧烈跳动的点。
5.2 奇异值拐点不明显导致阶次误判
现象:对块 Toeplitz 矩阵做 SVD 后,奇异值从第 1 阶到第 50 阶均匀下降,找不到明显台阶,无法直接断定“真实阶次是 8”之类。
原因:响应信噪比低,噪声奇异值和信号奇异值混在一起;块行数 i 取值太小,统计信息不足;采集时长太短,协方差长滞后项估值方差大。
解决:不要试图从奇异值曲线上定阶。把最大扫描阶次设高(比如 i·m/2 的 80%),依靠稳定图来定阶。同时增加块行数 i,比如从 20 提到 40,再看奇异值曲线——如果下降斜率在某个区间变平,那附近就是信号子空间边界。实际操作中,我更依赖“稳定点数量随阶次变化的曲线”:当阶次增长到真实阶次附近时,稳定点数会突增,继续增加阶次则稳定点数不再明显上升,这个“平台”就是最佳阶次范围。
5.3 采样率与频带不匹配,模态“消失”了
现象:稳定图目标频段里只有非常稀疏的稳定点,明明结构有 8 Hz 模态,图上却在 7~9 Hz 之间一片空白。
原因:采样率过高导致目标模态在总频率范围内占比太小,SVD 的主能量被低频或高频噪声抢走;采样率过低则模态频率接近奈奎斯特频率,离散状态矩阵产生混叠。
解决:先对数据做功率谱估计,确认目标模态频率相对采样率的位置。一般让目标最高模态频率小于采样率的 1/5,例如 100 Hz 采样率,最高关注模态 15 Hz 左右,余量足够。如果目标模态集中在低频,比如 0.2~2 Hz,但原始记录是 200 Hz 采样,应该先降采样到 20~30 Hz 再跑 SSI,否则 Toeplitz 里大量高频噪声成分会淹没低频模态。降采样前务必做抗混叠滤波。
5.4 传感器数量不足导致振型畸变
现象:稳定图上频率、阻尼都稳定,但振型形状不连续,MAC 匹配率极低,相邻测点的振型符号突变。
原因:测点数量少,空间采样率不足,振型发生空间混叠;参考点选在响应很小或噪声很大的位置,导致投影方向失真;局部模态没有对应测点,振型分量缺失。
解决:硬件条件允许时,至少布设 3~5 个测点覆盖结构主要变形形态;测点避开节线位置。计算振型前对模态向量做归一化,例如以幅值最大通道为 1,其他通道取比例。如果只能有少量测点,那么不要强求 MAC 判据,改用频率和阻尼双判据定模态,振型仅作为参考。传感器数 m 如果小于 3,稳定图的可靠性会显著下降,这时优先考虑频域单参考点方法或用多批次数据平均。
6. 把稳定图做成自动化流水线:一个可以直接落地的技巧
稳定图人工判读非常耗时间,但完全自动化又容易把噪声模态当真。我的落地习惯是“自动化粗筛 + 定向复核”。先跑一遍自动频率直方图候选,再对每个候选频率点检查两个指标:稳定点命中率和阻尼离散度。命中率越高、阻尼离散度越低,该候选越可信。
下面这段代码可以直接用在你的项目里。它接收稳定点数据,输出一个带评分列表的候选模态:
def rank_stable_candidates(freqs, stable_flags, min_hits=3, bin_width=0.02): import numpy as np freq_bins = np.arange(0, np.max(freqs)+bin_width, bin_width) hit_counts, edges = np.histogram(freqs[stable_flags], bins=freq_bins) candidates = edges[:-1][hit_counts >= min_hits] scored = [] for f in candidates: mask = (strong_stable_freqs := None) # 占位,实际按频率bin过滤 # 选中该bin附近0.1Hz内的所有稳定点阻尼 near = np.abs(freqs[stable_flags] - f) < 0.05 hit = near.sum() damp_vals = stable_damps[stable_flags][near] # 需要预先收集对应阻尼 # 阻尼离散度用标准差/均值,小于0.5认为是稳定 score = hit * (1 - np.std(damp_vals) / (np.mean(damp_vals)+1e-6)) scored.append((f, hit, score)) return sorted(scored, key=lambda x: -x[2])注意上面代码中strong_stable_freqs处是占位逻辑,实际项目里你需要同时保留每个稳定点的频率和阻尼,按频率分箱后计算命中率。评分公式里,“命中数×阻尼稳定度”兼顾了重复性和可靠性。
有了评分后,我会把排名前 N 个候选频率作为“备选模态”,再回到稳定图上,放大该频率附近 ±0.2 Hz 的窗口,人工确认是否存在成柱的稳定点。如果确认,就用该频率附近稳定模态的均值作为最终结果。这个流程把机器的批量能力和人的判断力结合起来,既不怕漏点,也不怕误报。
最后说一个我踩过多次的教训:SSI 稳定图不是一次跑完就能交差的,必须用不同块行数 i 和不同数据段做交叉验证。我习惯把同一组 120 秒数据切成长度 60 秒的两段,分别跑 SSI,若两段稳定图在同一频率出现稳定柱,才把这个模态认定为真实。否则,即便稳定图再漂亮,也可能是这段数据特有的噪声巧合。希望这个习惯能帮你少走弯路,也希望这篇拆解能让你在下次打开SSI.rar时,直接照着把稳定图跑出来,并从里面取到可信的模态参数。
本文还有配套的精品资源,点击获取