简介这份资源聚焦欠定盲源分离UBSS这一信号处理方向面向音频处理、医学成像、通信解码及金融数据分析等领域的学习者与研究人员帮助理解观测通道少于源信号数量时的分离思路。压缩包共5个文件以4个m脚本文件和1个txt说明文件为主脚本涵盖SCAN、DBSE、sobi、vecsymm等算法实现readme则提供使用说明整体约5KB体量轻便适合快速上手与二次修改。资源围绕独立成分分析、稀疏表示、联合诊断学习及机器学习等主流方法展开可配合音频、脑电等场景验证分离效果。目前已有254人学习下载适合作为算法入门与实验对照的参考材料便于读者在理解混合模型与线性系统的基础上进一步探索噪声抑制、非高斯假设限制等挑战提升盲信号分离的实践能力。1. SCAN 欠定盲源分离当观测通道比源信号还少时怎么把信号拆开做设备故障诊断或者无线频谱监测的工程师迟早会撞上这个场景手上只有两三个传感器现场却混着四五个甚至更多独立信号源想把它们一一拆出来。这就是欠定盲源分离要解决的问题——观测通道数少于源信号数传统独立成分分析直接失效因为混叠矩阵连方阵都不是求逆这条路根本走不通。SCAN 这类基于稀疏分量分析的方案核心思路是利用信号在某个变换域里的稀疏性先估计混叠矩阵再逐个恢复源信号。它适合做机械振动多故障源分离、脑电多通道伪迹剔除、通信信号盲解调前处理的从业者。下面把我实际落地时踩过的路讲清楚从数据模型到代码复现到参数调优能直接抄。2. 欠定盲源分离的数学底子与 SCAN 的切入点2.1 为什么观测通道少于源信号就变成了病态问题先把这个问题的数学形式摆出来。标准线性瞬时混合模型是x(t) A·s(t) n(t)其中 x(t) 是 M 维观测向量s(t) 是 N 维源向量A 是 M×N 的混叠矩阵。适定和超定情况下 M ≥ NA 是方阵或高矩阵可以用独立成分分析ICA通过最大化非高斯性来估计 A 的逆。但欠定情况下 M NA 是扁矩阵它的零空间非平凡意味着有无数个 s(t) 能产生同一组观测 x(t)。这就是欠定盲源分离的根本困难解不唯一。那为什么还能做关键在于引入额外假设。最常用也最有效的假设是稀疏性——源信号在某个变换域时域、频域、时频域中任意时刻只有少数几个源是活跃的。在这个假设下虽然全局解不唯一但每个时刻的局部混合可以看成是少数几个源列的线性组合几何上表现为观测点聚集在由 A 的列向量张成的低维子空间里。SCANSparse Component Analysis via clustering and Nullspace也有文献叫 Sparse Component Analysis Network正是利用这个几何结构先通过聚类估计混叠矩阵 A 的列方向再用最短路径或线性规划恢复源信号。2.2 SCAN 方法的两阶段框架混叠矩阵估计与源恢复SCAN 的流程分两大阶段每个阶段都有独立的算法选择和参数敏感性。第一阶段混叠矩阵估计。对观测信号做短时傅里叶变换STFT得到时频域表示 X(t,f)。在稀疏假设下每个时频点上只有一个源主导此时观测向量 X(t,f) 的方向就近似等于 A 中对应列向量的方向。把所有时频点的方向向量画在单位超球面上它们会聚成 N 个簇每个簇的中心方向就是 A 的一列。聚类算法常用 K-means 或者更鲁棒的 DBSCAN后者不需要预先知道源数目对噪声点也有更好的容忍度。第二阶段源信号恢复。有了 A 的估计值  之后每个时频点上求解一个欠定方程 x ·s。由于稀疏性可以用 L1 范数最小化来求最稀疏解min ||s||₁ subject to x ·s。这本质上是一个线性规划问题可以用基追踪Basis Pursuit或者匹配追踪Matching Pursuit系列算法求解。最后对每个源的时频表示做逆 STFT得到时域波形。整个框架听起来清晰但实际落地时每一步都有坑。下面用 Python 代码把最小可复现版本跑通。2.3 用 Python 跑通 SCAN 的最小复现流程先安装依赖核心是 numpy、scipy 和 scikit-learn用于聚类pip install numpy scipy scikit-learn matplotlib然后生成模拟数据并执行完整的 SCAN 流程import numpy as np from scipy.signal import stft, istft from sklearn.cluster import KMeans np.random.seed(42) # 1. 生成模拟源信号 fs 1000 # 采样率 Hz T 2.0 # 时长秒 t np.arange(0, T, 1/fs) N 4 # 源信号数欠定源多于观测 # 四个源两个正弦、一个方波、一个调频 s1 np.sin(2 * np.pi * 50 * t) s2 np.sin(2 * np.pi * 120 * t) s3 np.sign(np.sin(2 * np.pi * 30 * t)) s4 np.sin(2 * np.pi * (200 100 * t) * t) S np.array([s1, s2, s3, s4]) # shape: (4, len(t)) # 2. 构造欠定混合 M 2 # 观测通道数M N欠定 A_true np.array([[0.8, 0.3, 0.5, 0.2], [0.2, 0.7, 0.4, 0.9]]) # 2x4 混叠矩阵 X A_true S 0.01 * np.random.randn(M, len(t)) # 加少量噪声 # 3. STFT 变换到稀疏域 f, t_stft, Zxx stft(X, fsfs, nperseg128, noverlap96) # Zxx shape: (M, freq_bins, time_frames) # 4. 提取方向向量并聚类估计混叠矩阵 # 取能量足够大的时频点避免噪声主导方向 power np.sum(np.abs(Zxx)**2, axis0) # (freq_bins, time_frames) threshold np.percentile(power, 70) # 只保留能量前30%的点 mask power threshold # 提取方向向量归一化到单位球面 directions [] for fi in range(Zxx.shape[1]): for ti in range(Zxx.shape[2]): if mask[fi, ti]: v Zxx[:, fi, ti] norm np.linalg.norm(v) if norm 1e-10: directions.append(v / norm) directions np.array(directions) # shape: (num_points, M) # K-means 聚类簇数设为源数 N kmeans KMeans(n_clustersN, n_init10, random_state42) labels kmeans.fit_predict(directions) A_est kmeans.cluster_centers_.T # shape: (M, N) # 修正列顺序和符号与真实 A 对齐仅用于评估 from scipy.optimize import linear_sum_assignment cost np.abs(A_est.T A_true) row_ind, col_ind linear_sum_assignment(-cost) A_est A_est[:, col_ind] print(真实混叠矩阵 A:\n, A_true) print(估计混叠矩阵 Â:\n, np.round(A_est, 3))这段代码的关键逻辑说明STFT 参数选择nperseg128对应频率分辨率约 7.8 Hznoverlap96给出 75% 重叠保证时频面上有足够密的采样点用于聚类。窗长太短频率分辨率不够太长则时域稀疏性被破坏。能量阈值np.percentile(power, 70)只保留能量最高的 30% 时频点。这是欠定盲源分离里非常关键的一步——低能量点基本是噪声方向向量随机散布会严重干扰聚类。阈值太高则有效点太少聚类不稳定。方向向量归一化把每个时频点的观测向量投影到单位球面上这样聚类只关心方向不关心幅度符合稀疏假设下“单源主导”的几何直觉。K-means 的局限这里用了 K-means 是因为源数已知。实际场景中源数往往未知需要用 DBSCAN 或者基于特征值衰减的估计方法。另外 K-means 对初始值敏感n_init10是必要的。跑完这一步你会看到估计的混叠矩阵和真实值已经比较接近。接下来是源恢复。# 5. L1 范数最小化恢复源信号 from scipy.optimize import linprog def l1_min_recover(x_obs, A_est): 对单个时频点求解 min ||s||_1 s.t. A_est s x_obs M, N A_est.shape # linprog 标准形式: min c^T z, s.t. A_eq z b_eq, bounds # 变量拆分: s u - v, u,v 0, |s| u v c np.ones(2 * N) A_eq np.hstack([A_est, -A_est]) b_eq x_obs.real # 先处理实部 res linprog(c, A_eqA_eq, b_eqb_eq, bounds[(0, None)] * (2 * N), methodhighs) if res.success: s res.x[:N] - res.x[N:] return s return np.zeros(N) # 对每个时频点恢复 S_est_stft np.zeros((N, Zxx.shape[1], Zxx.shape[2]), dtypecomplex) for fi in range(Zxx.shape[1]): for ti in range(Zxx.shape[2]): x_point Zxx[:, fi, ti] s_real l1_min_recover(x_point.real, A_est) s_imag l1_min_recover(x_point.imag, A_est) S_est_stft[:, fi, ti] s_real 1j * s_imag # 6. 逆 STFT 恢复时域波形 S_est np.zeros((N, len(t))) for i in range(N): _, s_rec istft(S_est_stft[i], fsfs, nperseg128, noverlap96) S_est[i, :len(t)] s_rec[:len(t)] # 计算恢复信号与真实源的相关系数 for i in range(N): corr np.corrcoef(S_est[i], S[i])[0, 1] print(f源 {i1} 恢复相关系数: {corr:.4f})这段恢复代码有几个实操要点逐点线性规划非常慢。上面这个双重循环在 2 秒数据上可能要跑几分钟。实际工程中会用匹配追踪OMP或者迭代硬阈值IHT来加速牺牲一点精度换十倍以上的速度。实部和虚部分开处理是一种简化。更严谨的做法是把复数线性规划转成实数形式统一求解但代码复杂度会上升。methodhighs是 scipy 目前最快的线性规划求解器比默认的 simplex 快很多。如果 scipy 版本较老没有这个选项升级到 1.6 以上。相关系数是评估恢复质量的直接指标。如果某个源的相关系数低于 0.8通常意味着聚类阶段把两个源的方向混在了一起或者该源在时频域不够稀疏。3. 混叠矩阵估计的精度决定一切聚类参数怎么调3.1 源数目未知时用 DBSCAN 替代 K-means上一章用了 K-means前提是知道源有几个。实际项目中源数往往未知这时候 DBSCAN 是更务实的选择。它不需要预设簇数还能把噪声点标记为 -1 直接丢弃。代价是要调两个参数邻域半径eps和最小样本数min_samples。from sklearn.cluster import DBSCAN from sklearn.neighbors import NearestNeighbors # 用 k-距离图辅助确定 eps k 5 nbrs NearestNeighbors(n_neighborsk).fit(directions) distances, _ nbrs.kneighbors(directions) k_dist np.sort(distances[:, -1]) # 找 k-距离曲线的拐点作为 eps # 实际使用时画出 k_dist 曲线肉眼找拐点这里用差分近似 diff np.diff(k_dist) eps_auto k_dist[np.argmax(diff) 1] print(f自动估计的 eps: {eps_auto:.4f}) db DBSCAN(epseps_auto, min_samples10) labels_db db.fit_predict(directions) n_clusters len(set(labels_db)) - (1 if -1 in labels_db else 0) print(fDBSCAN 发现 {n_clusters} 个簇) # 用各簇均值作为混叠矩阵列向量 A_est_db np.zeros((M, n_clusters)) for i in range(n_clusters): cluster_points directions[labels_db i] center cluster_points.mean(axis0) A_est_db[:, i] center / np.linalg.norm(center)eps的选取是 DBSCAN 最玄学的地方。k-距离图是个靠谱的辅助工具把每个点到第 k 近邻的距离排序画出来曲线拐点对应的距离就是合适的eps。min_samples一般取观测维度的 2 到 3 倍比如 M2 时取 5 到 10。如果聚类结果里 -1 标签占比超过 30%说明eps太小大量有效点被当成噪声丢掉了。3.2 时频点筛选能量阈值和单源主导判据方向向量聚类的前提是“每个时频点只有一个源主导”。但现实中这个条件并不总是满足尤其在两个源频率重叠的区域。除了能量阈值还可以加一个单源主导判据计算该时频点的观测向量与最近簇中心的余弦相似度低于某个阈值比如 0.95就丢弃。# 在聚类完成后用簇中心反筛时频点 cos_threshold 0.95 valid_points [] valid_indices [] for idx, (d, label) in enumerate(zip(directions, labels)): if label -1: continue center kmeans.cluster_centers_[label] cos_sim np.dot(d, center) / (np.linalg.norm(d) * np.linalg.norm(center)) if cos_sim cos_threshold: valid_points.append(d) valid_indices.append(idx) print(f原始有效点: {len(directions)}, 筛选后: {len(valid_points)})这个二次筛选能显著提升混叠矩阵估计精度代价是可用于恢复的时频点变少可能影响源恢复的完整性。实践中cos_threshold取 0.9 到 0.98 之间根据源之间的频率重叠程度来定。重叠严重就放宽到 0.9源之间频率分得很开可以严到 0.98。3.3 混叠矩阵列顺序模糊与幅度歧义的工程处理盲源分离有两个固有的模糊性排列模糊和尺度模糊。排列模糊是指估计出来的 A 的列顺序和真实源顺序对不上尺度模糊是指每列可以乘以任意非零常数而观测不变。这两个模糊性在数学上无法消除但工程上可以处理。排列模糊的常见解法是如果你有每个源的先验频率范围或者波形特征可以在恢复后对每个源做频谱分析按特征匹配重新排序。尺度模糊则更简单——大多数应用只关心波形形状不关心绝对幅度把每个恢复源归一化到单位方差即可。如果确实需要绝对幅度可以用一个已知幅度的参考源做标定。# 按主频排序解决排列模糊 from scipy.fft import fft, fftfreq def dominant_freq(signal, fs): 返回信号的主频 spectrum np.abs(fft(signal)) freqs fftfreq(len(signal), 1/fs) positive freqs 0 return freqs[positive][np.argmax(spectrum[positive])] # 假设已知真实源的主频排序 true_freqs [dominant_freq(s, fs) for s in S] est_freqs [dominant_freq(S_est[i], fs) for i in range(N)] # 按主频匹配重排 order np.argsort(est_freqs) S_est_sorted S_est[order] print(重排后各源主频:, [f{dominant_freq(S_est_sorted[i], fs):.1f} for i in range(N)])这段代码展示了用主频特征解决排列模糊的思路。实际项目中特征可以更丰富谱峰位置、包络形状、循环平稳特征等。关键是找到对尺度变化不敏感、对源身份有区分度的特征。4. 源恢复阶段的加速与精度权衡4.1 匹配追踪替代线性规划快十倍但精度掉多少逐点线性规划在长数据上不可接受。匹配追踪Matching Pursuit系列是更实用的选择。正交匹配追踪OMP每次迭代选一个与残差最相关的原子然后正交投影更新残差。在欠定盲源分离里“原子”就是混叠矩阵的列向量。def omp_recover(x_obs, A_est, sparsity2): 正交匹配追踪恢复sparsity 为预期活跃源数 M, N A_est.shape residual x_obs.copy() indices [] s_est np.zeros(N, dtypecomplex) for _ in range(sparsity): # 找与残差最相关的列 correlations np.abs(A_est.conj().T residual) # 排除已选列 for idx in indices: correlations[idx] 0 new_idx np.argmax(correlations) indices.append(new_idx) # 最小二乘更新 A_sub A_est[:, indices] s_sub np.linalg.lstsq(A_sub, x_obs, rcondNone)[0] residual x_obs - A_sub s_sub for i, idx in enumerate(indices): s_est[idx] s_sub[i] return s_est # 用 OMP 替换线性规划速度对比 import time start time.time() S_omp np.zeros((N, Zxx.shape[1], Zxx.shape[2]), dtypecomplex) for fi in range(Zxx.shape[1]): for ti in range(Zxx.shape[2]): S_omp[:, fi, ti] omp_recover(Zxx[:, fi, ti], A_est, sparsity2) print(fOMP 耗时: {time.time() - start:.2f} 秒)sparsity参数控制每次恢复时假设有几个源同时活跃。在稀疏假设下这个值通常取 1 到 3。取太大会引入虚假源取太小会丢失同时活跃的源。一个经验做法是先用 1 跑一遍看恢复信号的相关系数如果某些源质量差就增加到 2 再试。OMP 的速度比线性规划快一个数量级以上但精度在源之间相关性高时会下降。如果对精度要求极高且数据量不大线性规划仍然是金标准。4.2 时频掩膜法更简单但依赖混叠矩阵精度如果混叠矩阵估计得足够准还有一个更快的路子时频掩膜。对每个时频点计算它属于每个源的概率基于到各簇中心的距离然后做软掩膜或者硬掩膜。def tf_mask_recover(Zxx, A_est): 基于时频掩膜的源恢复 M, N A_est.shape F, T Zxx.shape[1], Zxx.shape[2] S_masked np.zeros((N, F, T), dtypecomplex) for fi in range(F): for ti in range(T): x_point Zxx[:, fi, ti] # 计算到每个源方向的投影能量 energies np.zeros(N) for n in range(N): a_n A_est[:, n:n1] proj a_n np.linalg.pinv(a_n) x_point energies[n] np.linalg.norm(proj)**2 # 软掩膜按能量比例分配 total energies.sum() if total 1e-10: mask energies / total for n in range(N): S_masked[n, fi, ti] mask[n] * x_point np.linalg.pinv(A_est[:, n:n1]) A_est[:, n:n1] return S_masked时频掩膜法的恢复质量高度依赖混叠矩阵的精度。如果 A 估计有 5% 以上的误差掩膜就会把源之间的能量搞混恢复信号里出现明显的串扰。它的优势是速度快没有迭代求解适合做实时处理的预处理。4.3 恢复质量评估相关系数、信干比与主观听感评估恢复质量不能只看一个指标。相关系数反映波形形状的相似度但对幅度不敏感。信干比SIR反映源之间的串扰程度。如果是语音信号还得靠人耳听。def evaluate_separation(S_true, S_est): 计算各源的相关系数和信干比 N S_true.shape[0] results [] for i in range(N): # 相关系数 corr np.corrcoef(S_true[i], S_est[i])[0, 1] # 信干比目标源能量 / 其他源泄漏能量 # 用最小二乘拟合找泄漏系数 coeffs np.linalg.lstsq(S_true.T, S_est[i], rcondNone)[0] target_energy (coeffs[i] * S_true[i])**2 interference sum((coeffs[j] * S_true[j])**2 for j in range(N) if j ! i) sir 10 * np.log10(target_energy.sum() / (interference.sum() 1e-10)) results.append({source: i1, corr: corr, SIR_dB: sir}) return results for r in evaluate_separation(S, S_est): print(f源 {r[source]}: 相关系数{r[corr]:.4f}, SIR{r[SIR_dB]:.2f} dB)相关系数高于 0.9 且 SIR 高于 15 dB 通常算合格。如果 SIR 低于 10 dB说明串扰严重需要回头检查混叠矩阵估计或者增加稀疏性约束。5. 避坑与排查欠定盲源分离翻车实录5.1 聚类簇数设错导致混叠矩阵完全跑偏现象恢复出来的信号和真实源完全不相关相关系数在 0.2 以下但代码没有报错。原因K-means 的n_clusters设成了错误的源数。如果实际有 4 个源但设成 3K-means 会把两个方向相近的源合并成一个簇导致混叠矩阵少一列后续恢复时这两个源被当成一个处理。设成 5 则多出一个虚假簇恢复时引入不存在的源。解决源数未知时用 DBSCAN或者用特征值分解方法估计源数——对观测协方差矩阵做特征值分解看特征值衰减的拐点。另一个实用技巧是跑多个不同簇数的 K-means比较聚类轮廓系数silhouette score选最高的那个。5.2 时频域稀疏性不足导致方向向量弥散现象方向向量在单位球面上不是聚成几个清晰的簇而是弥散成一片聚类轮廓系数低于 0.3。原因源信号在所选变换域里不够稀疏。比如两个源在大部分时间频率上同时活跃每个时频点都是多源混合方向向量就失去了“单源主导”的几何意义。常见于源之间频率重叠严重或者源本身是宽带非稀疏信号。解决换变换域。STFT 不稀疏就试小波变换或者同步挤压变换。如果源在时域有稀疏性比如脉冲信号直接在时域做聚类。另一个方向是增加观测通道数——哪怕从 2 个增加到 3 个稀疏性要求都会大幅降低。5.3 噪声能量估计偏差让阈值筛选失效现象能量阈值筛选后有效时频点数量波动很大换一段数据聚类结果就完全不同。原因用了固定的百分位阈值比如 70%但不同数据段的噪声水平不一样。噪声大的段前 30% 的点里也混着大量噪声噪声小的段可能把有效信号也筛掉了。解决改用自适应阈值。先估计噪声底——取时频能量分布的最低 10% 分位作为噪声水平然后阈值设为噪声水平的 3 到 5 倍。这样阈值随数据自适应鲁棒性好很多。5.4 线性规划恢复时源信号被过度稀疏化现象恢复出来的源信号波形被“压扁”了峰值被削平听起来有失真。原因L1 范数最小化倾向于把小的系数压到零导致恢复信号过于稀疏。这在源信号本身不是严格稀疏时尤其明显。解决改用 L1-L2 混合范数或者加一个小的 L2 正则项。另一个办法是减小sparsity参数允许更多源同时活跃。如果对波形保真度要求高考虑用迭代硬阈值IHT替代线性规划它在稀疏度和保真度之间平衡得更好。5.5 逆 STFT 时相位不一致导致重建波形畸变现象恢复的时域波形在拼接处有跳变频谱上出现不该有的高频分量。原因对每个时频点独立做 L1 恢复时相邻时频点的相位没有连续性约束。逆 STFT 时相位不连续会导致重建波形出现伪影。解决在恢复阶段加相位平滑约束或者用重叠相加法时确保窗函数满足 COLAConstant Overlap-Add条件。更简单的做法是增大 STFT 重叠率到 87.5%noverlap 7/8 * nperseg让相位不连续的影响被平滑掉。6. 从仿真到实测SCAN 在机械振动信号上的调参心得仿真跑通只是第一步实测信号才是真正的考验。我拿一个双通道加速度传感器采的齿轮箱振动信号试过里面混着轴承故障冲击、齿轮啮合振动和电机电磁噪声三种源。仿真里表现很好的参数实测时几乎全要重调。第一个要改的是 STFT 窗长。仿真里 128 点窗对应 7.8 Hz 分辨率够用但实测轴承故障冲击的重复频率只有十几 Hz7.8 Hz 的分辨率根本分不开。我把窗长加到 512 点频率分辨率降到 2 Hz冲击的边频带才在时频面上清晰起来。代价是时域分辨率变差但振动信号本来就是频域特征更关键。第二个要改的是能量阈值。实测信号的噪声底比仿真高得多固定百分位阈值筛出来的点里混了大量噪声。我改成先对时频能量做对数变换然后用 Otsu 方法自动找双峰之间的谷底作为阈值。这个自适应方法在不同工况下都稳定不用每次手动调。第三个要改的是聚类距离度量。仿真里用欧氏距离就够了但实测信号的方向向量在球面上分布不均匀欧氏距离会把密集区域的点过度合并。换成余弦距离后簇的分离度明显改善。scikit-learn 的 KMeans 不支持自定义距离得用sklearn.cluster.SpectralClustering或者自己实现。最后一个血泪教训混叠矩阵的估计精度比源恢复算法重要得多。我花了大量时间调 L1 求解器的参数效果提升有限。后来回头把聚类环节的时频点筛选和距离度量优化了一下混叠矩阵的列方向误差从 8 度降到 2 度恢复信号的 SIR 直接涨了 6 dB。所以如果你的 SCAN 效果不好先检查混叠矩阵估计别急着换恢复算法。实测中还有一个仿真里遇不到的问题传感器通道之间的时间延迟。两个加速度计安装位置不同同一个冲击到达两个传感器有时间差这在频域表现为相位随频率线性变化。如果不补偿混叠矩阵的估计会有系统性偏差。补偿方法是在 STFT 后对每个通道乘以一个相位校正因子校正量通过互相关估计。这套流程跑通之后我把代码封装成了一个可配置的管道STFT 参数、聚类参数、恢复算法都做成配置文件。换一个工况只需要改配置不用动代码。现在处理一段 10 秒的双通道振动信号从 STFT 到恢复出三个源波形全程大约 15 秒其中聚类占 2 秒OMP 恢复占 12 秒。如果换成线性规划同样数据要跑 8 分钟以上。所以除非精度实在不够我一般默认用 OMP。希望帮到你。本文还有配套的精品资源点击获取