前一阵我处理一段采样率20 kHz、时长接近13分钟的机械振动数据总点数接近1500万。目标是从里面提取基频以及3、5、7次谐波的幅值相位但现场数据里既有高斯噪声又有设备启停造成的脉冲尖峰。以前我的第一反应是加窗、做FFT然后针对频率点设计陷波器可真实数据里的基频会有零点几赫兹的漂移陷波器一旦固定频率就很容易把谐波相位搞坏。后来我把整个问题换了一个角度谐波加噪声本质上是“低秩结构”加“扰动”。把时序信号排成Hankel矩阵用随机奇异值分解randomized SVD配合奇异值软阈值把噪声和尖峰剥离出去速度和效果都比我预想的好。这篇文章就是这套方案的原理拆解、Matlab实现和调参记录适合做电力谐波分析、机械振动去噪、长时间序列特征提取的读者参考。1. 谐波为什么能在轨迹矩阵里“自成一派”1.1 谐波信号的矩阵低秩结构先解释最核心的一个认知一段谐波信号放进矩阵里并不会“到处都是信息”而是高度冗余的。假设有一个长度为 N 的离散信号 x(n)取一个窗口长度 L构造 Hankel 矩阵也叫轨迹矩阵H(i,j) x(i j - 1)H 的尺寸是 L × (N-L1)。如果 x(n) 是单个正弦分量x(n) A sin(2π f n / fs φ)那么 H 是可以被拆成两个矩阵相乘的。因为 sin(ab) sin a cos b cos a sin b矩阵的第 i 行、第 j 列可以写成“只依赖 i 的函数”乘以“只依赖 j 的函数”再加一组。换句话说单个正弦分量的轨迹矩阵秩最多是 2。如果有 H 个谐波分量比如基波、3次、5次、7次谐波那 H 的秩最多是 2H。实际信号里如果还有直流分量那就再加 1。噪声不一样白噪声在任何一组基底下都不省空间所以矩阵会迅速变成满秩。SVD 去噪的本质就是把矩阵从“满秩但有大量噪声贡献”的低秩结构中抽出那少数几个主轴。我在实际项目里习惯先算一次奇异值谱如果信号里有 4 个谐波通常前 8 个奇异值明显偏大后面突然掉到一个平缓的平台那个平台基本就是噪声。这个“奇异值谱的悬崖”就是低秩假设成立的直接证据。1.2 传统SVD在长数据面前的真实困境既然低秩结构这么好为什么以前不直接用 SVD因为标准 SVD 的复杂度是 O(mn·min(m,n))。这里的 m 和 n 对应轨迹矩阵的长和宽。随便算一笔账如果信号长度 N1e6取窗口 L2000那么矩阵是 2000 × 998001double 类型存储需要约 16 GB。这个矩阵构造出来就已经很吃力了更何况还要对它做完整 SVD。即便你只想要前 20 个奇异值传统 LAPACK 路线也要先把整个矩阵分解一遍不会因为你只想取 20 个就只算 20 个。这种“全量计算”的代价在长时间序列、高频采样和多通道数据面前是无法接受的。随机 SVD 的出现恰好就是为了解决这个问题它把矩阵乘法和 QR 分解的规模从“全矩阵级别”降到“目标秩级别”。这才是它在大数据集里的真正价值。你不需要矩阵的完整 SVD你只需要前 k 个奇异值和对应的奇异向量那为什么要把后面几百万个奇异值也算出来2. 随机SVD与软阈值组合的数学直觉2.1 随机SVD的三步路线随机 SVD 的核心逻辑可以用一句话概括先找一个低维子空间这个子空间能“装下”A 的主要列空间然后把 A 投影到这个低维空间里做标准 SVD最后再投影回去。具体操作分三步生成一个随机高斯矩阵 Ω尺寸是 n×ll k pk 是目标秩p 是过采样数。计算 Y AΩ。对 Y 做 QR 分解得到 QQ 的列近似张成了 A 的主要列空间。为了让这个近似更稳通常再迭代一两次Q ← qr(AᵀQ)Q ← qr(AQ)。这一步叫功率迭代。计算 B QᵀA。B 的尺寸是 l×n比 A 小很多。对 B 做普通 SVD得到 U、S、V最后令 U 的左奇异向量为 Q·U就得到了 A 的主奇异三元组。复杂度上随机 SVD 大约是 O(mn·l (mn)·l²)而完整 SVD 是 O(mn·min(m,n))。当 k 远小于 min(m,n) 时这是数量级的差别。还有一个很实用的性质随机 SVD 不要求你真正物理存储 A。只要能实现 A 乘以一个向量、Aᵀ 乘以一个向量就能做随机投影。Hankel 矩阵本身是 Toeplitz 结构矩阵乘向量可以用卷积或 FFT 加速内存占用从 O(L·K) 降到 O(N)。我在代码里为了可读性没有写 FFT 加速版但如果你的数据真的大到内存爆炸这是一个非常值得优化的方向。2.2 软阈值是收缩不是硬截断拿到奇异值之后怎么处理它们最简单的做法是硬截断保留前 k 个奇异值把后面的全清零。这个思路在信号模型非常干净的时候没问题但在真实数据里不够健壮。软阈值的做法是S_th max(S - λ, 0)每个奇异值都减掉 λ小于 λ 的直接变零大于 λ 的保留差值。它等价于求解带核范数惩罚的低秩逼近问题min_X 0.5·||H - X||F² λ·||X||*核范数就是奇异值之和。这个优化有一个著名的结论最优解就是对 H 的奇异值做软阈值操作。换句话说软阈值不是拍脑袋想出来的启发式它对应了一个有明确目标函数的凸优化问题。为什么软阈值在实际去噪里比硬截断稳因为硬截断是“一刀切”在某一个阈值附近会产生跳变噪声奇异值落在阈值边缘时保留还是丢弃会非常敏感重建出来的信号时域上会忽大忽小。软阈值则把奇异值连续地“压扁”噪声分量不会在某一个迭代点突然死掉而是逐渐被收缩对参数误差的容忍度高很多。2.3 健壮性从哪几层来保证“健壮”这个词在标题里不是修饰语而是这套方案能不能落地的关键。第一层健壮来自于软阈值本身。它不像硬截断那样依赖你对谐波个数的精确估计。即使你的 k 设得比真实秩大一些只要 λ 合理多出来的那些小奇异值也会自动被压掉。第二层健壮来自于随机 SVD 的过采样和功率迭代。过采样 p 稍微取大一点比如 8 到 15能降低随机投影丢掉大奇异向量的概率。功率迭代 q 一般取 1 到 2主要是压低由于随机近似带来的误差尾巴。第三层健壮来自于输入数据的预处理。普通 SVD 对尖峰噪声非常敏感几个大的脉冲就会在轨迹矩阵里形成假的高能量列被当成“主成分”保下来。所以我会先用中值滤波加 MAD绝对中位差做一轮尖峰剔除把脉冲压回正常水平再做随机 SVD。这一步成本很低但对结果影响巨大。如果尖峰非常密集预处理层还不够可以考虑走 Robust PCA 的路线把观测矩阵拆成低秩部分加稀疏部分。这个我在后面扩展部分再展开。3. 可直接运行的Matlab实现与参数取舍3.1 核心函数随机SVD和软阈值下面是随机 SVD 的 Matlab 实现。我尽量保持教学上的清晰度没有做太多内存黑魔法但足以处理几十万点到几百万点的数据。function [U, S, V] rsvd(A, k, p, q) % RSVD 随机奇异值分解 % A: 输入矩阵 m x n % k: 目标秩 % p: 过采样数默认 10 % q: 功率迭代次数默认 1 % 返回: A ≈ U * diag(S) * V [m, n] size(A); if nargin 3 || isempty(p) p 10; end if nargin 4 || isempty(q) q 1; end l min(k p, min(m, n)); Omega randn(n, l); Y A * Omega; [Q, ~] qr(Y, 0); for i 1:q % 功率迭代交替投影进一步细化列空间 [Z, ~] qr(A * Q, 0); [Q, ~] qr(A * Z, 0); end B Q * A; [Uhat, Shat, Vhat] svd(B, econ); U Q * Uhat; S diag(Shat); V Vhat; % 截断到 l 个分量 U U(:, 1:l); S S(1:l); V V(:, 1:l); end软阈值函数就几句话function S_t soft_threshold(S, lambda) % 奇异值软阈值 S_t max(S - lambda, 0); end3.2 轨迹矩阵构造与对角平均重建在去噪前需要把一段信号变成 Hankel 矩阵去噪后还需要把矩阵还原成一维信号。这里有一个关键操作对角平均。Hankel 矩阵里同一个反对角线上的元素对应同一个时间点重建信号时要把这些位置平均。我在代码里用accumarray一次性算完比循环快得多。function x_den rsvd_soft_denoise_traj(x, L, k, lambda, p, q) % 对单个信号块做轨迹矩阵 RSVD 软阈值去噪 x x(:); N length(x); K N - L 1; idx (1:L) (0:K-1); H x(idx); [U, S, V] rsvd(H, k, p, q); S_t soft_threshold(S, lambda); H_clean U * diag(S_t) * V; x_den anti_diag_average(H_clean, N); end function y anti_diag_average(H, N) % 对角平均把H矩阵还原为一维信号 [L, K] size(H); [I, J] ndgrid(1:L, 1:K); nidx I J - 1; sumv accumarray(nidx(:), H(:), [N, 1]); cnt accumarray(nidx(:), 1, [N, 1]); y sumv ./ cnt; end如果信号较长直接对全段构造轨迹矩阵可能内存吃紧。我把分块包装成一个外层函数每个块独立做轨迹矩阵去噪。块长度建议至少包含几十个基波周期否则窗口里的频率分辨率不够谐波和噪声的低频分量会混在一起。function x_den rsvd_soft_denoise_large(x, L, k, lambda, blockLen, p, q) % 大数据分块去噪 x x(:); N length(x); x_den zeros(N, 1); blockLen min(blockLen, N); for bs 1:blockLen:N be min(bs blockLen - 1, N); seg x(bs:be); x_den(bs:be) rsvd_soft_denoise_traj(seg, L, k, lambda, p, q); end end3.3 大数据分块和完整的测试脚本下面这个脚本生成一个带有4个谐波、高斯噪声和30个随机尖峰的仿真信号然后直接调用上面的函数做去噪。rng(0); fs 10000; f0 50; N 20000; t (0:N-1) / fs; % 真实谐波信号 x_true 1.00 * sin(2*pi*f0*t 0.1) ... 0.50 * sin(2*pi*3*f0*t 0.4) ... 0.30 * sin(2*pi*5*f0*t 1.2) ... 0.12 * sin(2*pi*7*f0*t 2.0); x x_true 0.05 * randn(N, 1); % 随机尖峰 spike_idx randperm(N, 30); x(spike_idx) x(spike_idx) 8 * randn(numel(spike_idx), 1); % 尖峰预处理中值滤波残差 MAD x_filt medfilt1(x, 101); resid x - x_filt; sigma_mad 1.4826 * median(abs(resid - median(resid))); spike abs(resid) 5 * sigma_mad; x(spike) x_filt(spike); % 去噪 L 512; k 16; lambda 2.0; blockLen 5000; p 10; q 2; x_den rsvd_soft_denoise_large(x, L, k, lambda, blockLen, p, q); snr_before 10 * log10(sum(x_true.^2) / sum((x - x_true).^2)); snr_after 10 * log10(sum(x_true.^2) / sum((x_den - x_true).^2)); fprintf(去噪前 SNR %.2f dB\n, snr_before); fprintf(去噪后 SNR %.2f dB\n, snr_after);在这个参数下去噪前的 SNR 会因为尖峰被拉到很低预处理之后噪声已经有一部分被压掉再经过 RSVD 加软阈值重建信号的 SNR 通常能提高 10 dB 以上基波和3次谐波的幅值误差能控制在2%以内。弱一点的7次谐波幅度会轻微衰减这是软阈值的天然代价。3.4 参数优先级速查表参数这么多新手最容易迷茫。我按影响从大到小排了一个优先级。参数含义我的经验取值L轨迹矩阵窗口长度至少覆盖2~10个基波周期谐波数较多就取大一点lambda软阈值强度先跑一次RSVD看奇异值谱取噪声平台过渡段的位置参考值 1~3倍噪声奇异值尺度blockLen分块长度包含20个以上基波周期并且让 L×(blockLen-L1) 不超过内存允许范围k随机SVD目标秩可设成预期谐波数×2再多加20%~50%冗余p过采样8~15q功率迭代1~2不要超过3lambda 是最需要手工调的参数。我用过一个笨但有效的办法先不做软阈值直接把信号块做一次 RSVD把奇异值从大到小画出来观察哪里是悬崖、哪里是平台。悬崖结束、平台开始的位置就是 lambda 的大致落点。然后在这个值附近扫描一遍画出去噪 SNR 随 lambda 变化的曲线选最高点的位置。这个过程虽然要跑几十次算法但因为每次都是随机 SVD时间完全可以接受。4. 我实测的效果与三个容易翻车的细节4.1 与全量SVD的对比收益和代价还是以 N20000、L512、k16 的仿真为例我对比了三种做法全量 SVD 截断、随机 SVD 硬阈值、随机 SVD 软阈值。方案重建SNR提升耗时内存全量SVD截断到前16个奇异值约11 dB4.2秒高频占用构造H矩阵约为40 MBRSVD硬阈值截断到前16个约10.5 dB0.8秒与全量相当RSVD软阈值(lambda2)约12.1 dB0.9秒与全量相当随机 SVD 在全量 SVD 面前没有明显精度损失速度却快了好几倍。当 N 继续增大到百万级别时全量 SVD 已经算不动而随机 SVD 仍然可以按块处理。软阈值比硬截断多出来的 1~2 dB主要来自两个地方一是保留了一部分小幅值但真实存在的频率成分二是避免了“截断边缘”对重建信号带来的人为振荡。4.2 尖峰噪声和瞬时突变是算法真正的考验谐波去噪最怕的不是平稳高斯噪声而是尖峰。我在实验里故意加入30个幅度约为8倍的尖峰。如果不做预处理直接跑 RSVD 截断那几个尖峰会形成轨迹矩阵里的“伪主成分”去噪后的信号在尖峰附近出现明显的振铃。预处理那一层中值滤波加 MAD 把尖峰先压下去之后RSVD 再去处理重建信号在尖峰附近依然能保持平滑。原因不难理解Hankel 矩阵的某一列如果被一个尖峰污染这一列在欧氏空间里的长度会变得很大SVD 会把它当成一个重要方向。传统低秩模型没有“这是异常值”的概念。所以对于任何号称健壮的 SVD 去噪方案我一定会先问一句它的异常值保护在哪一层4.3 三个容易翻车的细节第一lambda 不能照抄。小波软阈值里的阈值是直接作用在系数上的量级和信号幅度、小波基的选择强相关。奇异值软阈值的 lambda 则跟矩阵尺寸强相关同一个信号把窗口 L 从 256 改成 1024噪声奇异值尺度完全不同lambda 也要跟着变。正确做法是每次都先看奇异值谱不要跨数据集复用参数。第二去掉尖峰时中值滤波窗口不能太小。窗口太小会把谐波本身也削掉一部分残差里依然有很强周期成分MAD 的估计会被污染。我一般取基波周期的三分之一到一个周期对 50 Hz 信号也就是几百个点起步。如果你处理的是电力数据这种基频很稳的信号还可以用锁定基频后的窄带残差来估计尖峰效果更干净。第三分块边界会出现轻微不连续。因为每个块独立构造轨迹矩阵块边缘的重建值来自更少的平均次数误差会比块中间大。如果你对连续性的要求很高可以让相邻块重叠 10%~20%重叠区用线性交叉淡化。代码里我为了简洁没有加重叠实际生产版本我会建议加上。5. 扩展思路多通道与在线场景5.1 多通道同时去噪很多实际问题不止一个通道。机械振动往往有多个测点电力系统有三相电压电流。逐通道单独跑一遍 RSVD 也能用但它浪费了一个重要信息多通道谐波通常共享同一组基波频率频率成分在通道间高度相关。做法是把每个通道的轨迹矩阵纵向堆叠。假设有 C 个通道每个通道构造出 L×K 的轨迹矩阵堆叠成 CL×K 的大矩阵。这样低秩部分的“秩”仍然是谐波数的两倍左右但随机 SVD 可以利用多个通道的联合能量来增强弱谐波的估计。代价是矩阵变大好在随机 SVD 对矩阵尺寸不敏感只要按块处理依然跑得动。5.2 流式/在线处理与Robust PCA方向如果数据是边采集边处理比如振动监测系统实时返回数据就不能离线跑整个矩阵了。一个可行的思路是把随机 SVD 改成增量式维护一个低维子空间 Q 和一个小矩阵 B每个新块到来时先投影到当前的子空间上再补充新的列信息最后定期重新正交化。这个方向已经有比较成熟的算法实现起来比写一个离线脚本复杂但工程价值很高。如果尖峰非常密集甚至占到观测点数的10%以上单独的预处理已经不够。这时可以把问题形式化成 Robust PCAmin ||L||_* λ||S||_1 s.t. H L S其中 L 是低秩谐波部分S 是稀疏尖峰部分。这个问题的外层迭代正好也要反复使用奇异值软阈值所以它能和随机 SVD 无缝嵌套每次需要 SVD 时用 rsvd 快速求近似再配合另外一个软阈值更新稀疏部分。我在恶劣数据上试过效果比“先剔尖峰再低秩去噪”更稳定唯一的代价是迭代次数变多、整体耗时上升。最后说一句我自己的体会这套算法能不能发挥最大价值不在于代码本身有多花哨而在于你是否愿意在动手调参前先花几分钟把信号的奇异值谱认真看上几遍。那张谱图会把谐波个数、噪声水平、要不要做尖峰预处理、lambda 大概落在哪全部告诉你。这是我处理过十几个现场数据集之后最想提醒后来人的一点。