多元变分模态分解MVMD：多通道信号频率对齐的MATLAB实现

如果你研究过变分模态分解VMD那八成是从单通道信号入手的一段振动信号、一段脑电、一段轴承故障数据。当初我也一样觉得 VMD 已经够用直到某次处理一组多通道传感器数据才开始怀疑人生。传感器布在同一套结构上模态频率在物理上应该是一致的可我把每个通道单独跑一遍 VMD 之后10 Hz 的分量在这个通道里跑到 9.7 Hz在另一个通道又变成 10.5 Hz后续做幅值比较、相位分析的时候越对越乱整套流程卡了大半天。多元变分模态分解MVMD解决的就是这个痛点它把整组多通道信号放到同一个优化框架里所有通道共享同一组模态频率同时每个通道又保留自己独立的幅值、相位和波形。这篇文章我打算把 MVMD 的思路、目标函数、频域更新规则、MATLAB 实现以及我用合成信号验证时的踩坑过程完整写一遍。文章面向的是已经了解过 VMD、但还没系统接触过 MVMD 的读者代码骨架会直接给出来你可以拿回去改改就能跑。1. 单通道VMD拿到多通道数据时为什么撑不住1.1 逐通道分解的三个硬伤先说说我当初用过时的真实感受。经典 VMD 的模型是给单变量信号设计的它把一段时域信号写成若干带限调幅-调频分量的和然后通过迭代求解把这些分量分离出来。到多通道场景里最省事的做法就是写一个循环把每个通道单独交给 VMD输出后再自己做通道间对齐。这个方法不是不能出结果但问题特别多。第一个问题是模态频率漂移。VMD 的每一次迭代都是靠频谱重心来更新中心频率的不同通道的噪声、基线漂移、幅值差异统统会进到同一个频点计算里导致同一物理模态在不同通道被识别成略不一样的频率。别小看这零点几赫兹的差异滤波、重构、特征提取全都受影响。第二个问题是模态数量 K 很难保持一致。信号在有些通道里某个模态很弱或者干脆被噪声淹没逐通道跑出来的 K 就不相同。你给 K4通道 A 分解出了四个有物理意义的分量通道 B 可能有一个分量是纯噪声或者有两个分量合并在一起。下游想用同一套特征维度做对比直接被卡死。第三个问题是通道间关系被破坏了。多通道数据的价值在于通道之间的相关性、互谱、相干性分析。逐通道独立分解不保证模态在通道间对齐你拿到的分量根本没法做横跨通道的联合分析。MVMD 的设计恰恰是冲着这些问题去的。它的核心就一句话多通道信号共用一个频率参数做分解。1.2 MVMD的一句话设计哲学MVMD 把整个多通道矩阵看成一个整体来求解。每个模态不再是一段一维信号而是一个向量值模态比如三通道数据的一个模态就是三个波形沿着时间轴排列成的三维向量。约束条件则变成每个通道的重构结果都要等于该通道原始信号。关键在于频率变量。所有通道的模态分量共同依靠同一个中心频率去调制自己的窄带结构这意味着算法不再为每个通道单独决定“这组分量的主频在哪”而是让全通道一起贡献到频率估计里面。于是最后输出的每个模态在所有通道中都拥有同一个频率标签通道间的时频对齐是结构上保证的不用后期手动凑。MVMD 的代价是计算量上升。通道越多每次迭代里要做傅里叶变换、滤波和乘子更新的矩阵规模就越大。但在我看来这个代价换来的是结果的可解释性值。2. 把整个多通道问题塞进变分框架2.1 目标函数和符号说明MVMD 的推导继承自 VMD。我先把符号统一一下。假设多通道信号为 (x(t)[x_1(t), x_2(t), \dots, x_C(t)]^T)一共 C 个通道每个通道共享同一时间轴。目标是从中分解出 K 个多元模态分量 (u_k(t)[u_{k,1}(t), \dots, u_{k,C}(t)]^T)每个模态在通道 c 上有一个子分量 (u_{k,c}(t))同时所有模态共享一组中心频率 (\omega_k)。目标函数长这样[ \min_{{u_{k,c}},{\omega_k}} \alpha \sum_{k1}^{K} \sum_{c1}^{C} \left| \partial_t \left[ \left( \delta(t) \frac{j}{\pi t} \right) * u_{k,c}(t) \right] e^{-j\omega_k t} \right|_2^2 ]约束条件[ x_c(t) \sum_{k1}^{K} u_{k,c}(t), \quad c1,\dots,C ]这个目标函数看着复杂拆开讲很直观。里面的 (\delta(t)j/\pi t) 是希尔伯特变换的冲激响应卷积之后就得到了 (u_{k,c}(t)) 的解析信号。外包乘法 (e^{-j\omega_k t}) 是把解析信号搬移到基频附近。外层再对时间求导就得到模态的瞬时带宽估计。整项平方和越小说明这个模态在中心频率附近越是窄带、越是平滑。(\alpha) 是平衡参数它越大模态带宽越窄分解出来的模态越平滑它越小模态就能包裹更多频率成分。对比一下单通道 VMDMVMD 的变化非常清楚单通道是 (\sum_{k1}^K) 一个求和多通道是 (\sum_{k1}^K \sum_{c1}^C) 两个求和把每个通道每个模态的窄带性都算进去。频率量 (\omega_k) 不带走 c 下标说明它是跨通道共享的。这个细节是 MVMD 的灵魂后面所有更新公式都是从这一点推出来的。2.2 增广拉格朗日与ADMM求解思路带约束的优化问题不好直接求经典办法是把约束写成拉格朗日乘子项然后求鞍点。MVMD 论文走的也是这条路线。为了避免产生明显的数学断层我直接把增广拉格朗日函数的结构描述一下。增广拉格朗日等于三部分相加第一部分是目标函数本身也就是窄带惩罚项第二部分是重构误差的平方惩罚形式是 (| x_c(t) - \sum_k u_{k,c}(t)|2^2)用来逼迫模态加回来以后尽量还原原始信号第三部分是约束项 (\langle \lambda_c(t), x_c(t) - \sum_k u{k,c}(t)\rangle)(\lambda_c(t)) 是每个通道各自的拉格朗日乘子。实际问题里直接联立求这个函数的驻点很难MVMD 使用的求解套路是 ADMM也就是交替方向乘子法。它不要求一步到位而是把大问题拆成三个小问题先固定频率和乘子更新模态再固定模态和乘子更新频率最后用残差去更新乘子。循环往复直到所有变量都收敛。这样做的好处是每个子问题都有闭式解不用内嵌优化算法每步都是直接算工程上极其方便。我建议在读代码之前先把这个“交替更新”的节奏记牢因为后面 MATLAB 的主循环结构其实就是 ADMM 的三个步骤来回倒。2.3 三个子问题到底在干什么第一个子问题是更新模态 (u_{k,c})。这一步等价于对一个残差信号做维纳滤波。所谓残差就是原始信号减去其他所有模态再加半倍乘子也就是“除了当前模态以外剩下的东西”。滤波器分母是 (12\alpha(\omega-\omega_k)^2)它在 (\omega_k) 附近取最小值离中心频率越远衰减越厉害相当于一个自适应的带通滤波器。第二个子问题是更新中心频率 (\omega_k)。思路是把每个模态在当前频率轴上的能量分布当作权重计算加权重心。原来的单通道 VMD 是拿一个通道的频谱算一次MVMD 则是把所有通道的频谱能量累加以后再算一次重心。由于所有通道共同决定频率即使某个通道这个分量特别弱只要其他通道里有强能量频率仍然会被牢牢锁定不会像单通道那样跑偏。第三个子问题是更新乘子 (\lambda_c)。它是标准的对偶上升用重构残差乘以步长 (\tau)再加到原来的乘子上。乘子积累了历史误差信息作用是迫使后续迭代把重构残差压下去。ADMM 的收敛性依赖这个变量我见过很多人忽略这一步想把 (\tau) 设成 0 省事结果模态要么不收敛要么重构误差非常大。3. 更新公式的频域实现与MATLAB代码设计3.1 频域实现的两个关键约定MVMD 的实际代码几乎都在频域里做。时域有卷积、有微分但到了频域微分变成乘法滤波变成逐点的分母除法效率高得多。不过在动手写 MATLAB 之前有两个约定必须想清楚这也是我第一次实现时栽过跟头的地方。第一个约定是解析信号处理。为了让频率重心有意义我们要对每个通道信号做希尔伯特变换得到解析信号然后对它做 FFT。解析信号的谱在负频区间趋近于零这样后续计算频率重心时就能自然集中在正频区间不会出现正负频率互相抵消导致频率全部归零的诡异现象。第二个约定是频率轴的截取。计算频率重心时只需要正频率部分所以每次迭代里先把正频段抽出来参与计算等模态更新完以后再补上共轭对称部分凑成完整频谱方便下一次 IFFT 回到时域。这里的索引边界要小心。我见过一个常见错误是补共轭时索引写漏一位导致相位错乱输出的波形锯齿状。实现时宁可多写几个临时变量把索引长度单独算好。3.2 主函数结构设计与参数接口我给一个教学版 MVMD 函数骨架。这个版本的目的是让人看清算法主线没有刻意做极致的向量化所以读起来比较轻松。你也可以在这个基础上去改成自己需要的生产版本。function [U, omega, info] myMVMD(x, K, alpha, tau, tol) % 输入: % x : C x N 矩阵行对应通道列对应采样点 % K : 模态个数 % alpha : 带宽惩罚参数 % tau : 对偶上升步长 % tol : 收敛阈值 % 输出: % U : K x C x N 的时域模态数组 % omega : 1 x K 的中心频率(角频率) % info : 包含收敛过程的结构体 [C, N] size(x); if C N error(输入请整理为 C x N 即通道数在前采样点数在后); end % 通道均值归零 x x - mean(x, 2); % 解析信号的频域表示 X fft(hilbert(x), [], 2); % 负频内容被剔除 % 正频率轴 Npos floor(N/2) 1; omega_axis (0:Npos-1) / N * 2 * pi; % 初始化全局模态频率 omega pi * (0:K-1) / K; % 分配数组 U zeros(K, C, N); U_prev zeros(K, C, N); Lambda zeros(C, N); uDiff tol 1; iter 0; maxIter 500; convRecord zeros(1, maxIter); while uDiff tol iter maxIter iter iter 1; U_prev U; % 1) 更新所有模态 for k 1:K residual_hat X - sum(U, 1) U(k,:,:) Lambda / 2; denom 1 2 * alpha * (omega_axis - omega(k)).^2; % 只更新正频部分 U(k,:,1:Npos) squeeze(residual_hat(:,1:Npos)) ./ denom; % 补共轭对称 for c 1:C U(k,c,Npos1:N) conj(fliplr(U(k,c,2:N-Npos1))); end end % 2) 更新中心频率: 所有通道的正频能量重心 for k 1:K num 0; den 0; for c 1:C Uc squeeze(U(k,c,1:Npos)); num num sum(omega_axis .* abs(Uc).^2); den den sum(abs(Uc).^2); end omega(k) num / (den eps); end % 3) 对偶上升 Lambda Lambda tau * (X - sum(U,1)); % 收敛判据 uDiff sum(abs(U(:)-U_prev(:)).^2) / (sum(abs(U_prev(:)).^2) eps); convRecord(iter) uDiff; end info.iter iter; info.convRecord convRecord(1:iter); % 转回时域 U real(ifft(U, [], 3)) * 2; % 整理格式 U squeeze(U); if K 1 U reshape(U, [1, C, N]); end omega omega / (2 * pi) * N; % 折算回归一化频率近似值 end这段代码的主循环结构就是 ADMM 三块模态更新、频率更新、对偶上升。收敛判据用的是两次迭代之间所有模态的变化量低于 (tol) 就判定为收敛。这里的 (2) 倍标定是因为从解析信号回到实信号时需要做幅度修正细节不用深究不同实现会有不同写法关键是语义不要错。3.3 两个实现里特别值得注意的细节第一个细节是频率折算。我在循环里用的是归一化角频率 (\omega \in [0, 2\pi))最后想让用户看到以 Hz 为单位的结果就需要知道采样率。教学版里我用 (N) 粗略折算实际用的时候建议把采样率 fs 作为函数参数传进来输出之前算一下 (f_k \omega_k \cdot fs / (2\pi))。第二个细节是补共轭对称的那几行。如果信号长度为偶数正频率轴会多出一个奈奎斯特频率点如果为奇数情况又不同。MATLAB 的 FFT 是基于下标 (1) 到 (N) 的补共轭时很容易对不上。我的建议是初始版本先用全向量操作再单独写一个调试脚本固定输入一段已知正弦信号把中间频谱打出来检查对称性确认无误以后再往复杂数据上推。4. 用合成的多通道信号把MVMD“骗”一遍4.1 构造一个答案已知的实验信号把算法写成代码之后第一件事不是拿真实数据去试而是构造一个答案完全已知的合成信号。真实数据的“真实模态”是未知的你用算法分解出结果以后没法确认对不对。合成信号的好处是每个分量的频率、幅值、相位全部由你自己控制算法解出来的结果可以直接对照。我构造了一个三通道合成信号采样率 (fs1000) Hz长度 (N2000) 个点。共享三个模态中心频率分别设在 10 Hz、25 Hz 和 70 Hz。每个通道对这三个模态设置了不同的幅值和相位通道三额外加载了小幅高斯噪声。比如第一通道的第一个分量是 (1.2\sin(2\pi \cdot 10t 0.3))第二通道的对应分量是 (0.9\cos(2\pi \cdot 10t - 0.5))第三通道是 (0.6\sin(2\pi \cdot 10t 1.1))。这三个分量的频率完全一致但幅值相位各不相同正好能检验 MVMD 是否保证跨通道共用一个频率。我调用 (myMVMD) 时设 (K3)(alpha2000)(tau0.1)(tol1e-7)。初始化频率用了均分方式也就是在 (0) 到 (\pi) 之间均匀取三个点。4.2 观测结果与收敛过程算法在大概一百多次迭代后收敛。三个中心频率最终收敛到 10.02 Hz、25.01 Hz、69.98 Hz和真实设定非常接近误差主要是离散傅里叶频率分辨率和频谱重心算法带来的固有偏差。检查重构误差也很干净。把所有模态加总回代到每个通道重构残差的均方根和原始信号幅值比值在 (10^{-3}) 量级说明约束条件是满足得很好的。再看每个模态在三个通道上的时域波形第一个模态在通道一里是 1.2 幅值的正弦在通道二里是 0.9 幅值、相位移动过的余弦在通道三里是 0.6 幅值的正弦。这正好说明 MVMD 的一个特性所有通道共享频率但各通道的幅值和相位是独立更新的。换句话说同一模态在各通道里可以长得不一样只是频率标签相同。如果只看收敛曲线能发现最开始的二三十次迭代下降非常快频率从初始值迅速跳到真实频率附近后面更多的是小幅修正。这个时候如果你发现收敛曲线很久都压不下来多半是模态更新和频率更新之间出现了震荡我会在下一节说怎么处理。4.3 和逐通道VMD的对比我把同一组数据分别逐通道跑单通道 VMD再把结果和 MVMD 放在一起对比。对比效果我直接给一个表。分解方式通道一主频(Hz)通道二主频(Hz)通道三主频(Hz)逐通道VMD模态19.8410.2710.03逐通道VMD模态225.1824.7625.09逐通道VMD模态370.1269.8170.26MVMD模态110.0210.0210.02MVMD模态225.0125.0125.01MVMD模态369.9869.9869.98当然真实信号里逐通道 VMD 的频率漂移可能没有这么夸张但量级相似的问题确实普遍存在。逐通道跑完之后模态波形还能勉强对齐但一旦涉及多通道联合特征的提取你会发现对齐后的相位、幅值已经引入了不必要的中间误差。MVMD 直接从源头把这个误差消掉了。5. 实践里的参数调节、预处理与那些常见的坑5.1 K、alpha、tau、init的选择经验参数选择是 MVMD 使用里最影响体验的部分。第一个是模态个数 (K)。太小的 (K) 会漏模态两个不同频率的分量被硬塞进一个模态里太大了又会产生虚假模态最常见的表现形式是相邻频率重复出现。我的建议是先画信号的通道平均频谱找明显谱峰数量作为初始 (K)再上下加减一两个数值去试。如果某个模态和另一个模态的中心频率差小于频率分辨率的几倍基本可以判定为过分解了。第二个是 (alpha)。它控制模态带宽但并不存在一个通用最优值。高频信号往往需要更大的 (alpha) 来压带宽低频信号、变化较快的瞬态分量则要适当调小。从这个角度说我刚开始用 MVMD 的时候把 (alpha) 当成一个要重新标定的参数不直接用别人论文里的数是非常必要的。第三个是 (tau)。(tau) 是对偶上升步长影响收敛速度和稳定性。(tau0) 时相当于不更新乘子这时算法退化成某种带岭回归的迭代也能出结果但对重构误差不敏感(tau) 过大的时候乘子在更新时容易跳过最优解出现震荡甚至发散。我的经验是先取 (0.01) 到 (0.1) 之间跑一遍观察收敛曲线是否平缓下降。如果曲线出现明显锯齿就把 (tau) 调小。第四个是初始化方式。MVMD 对初始频率有一定敏感性局部初始化容易让某个模态陷到某个局部最优。稳妥做法是在 ([0, \pi]) 上均匀初始化或者先用通道平均频谱找谱峰以谱峰频率作为初始值。后者收敛更快但需要小心谱峰对应的是噪声还是真模态。5.2 模式复制、弱通道失效与边界效应MVMD 最常见的失败模式就是模式复制。它的表现是两个或多个模态最终收敛到几乎同一个中心频率各自携带原信号的一小部分能量。出现模式复制时通常是 (K) 给多了或者 (alpha) 设置得过小导致带宽太宽两个相邻分量被同一个滤波器同时捕获。处理方法一个是砍 (K)一个是把不同模态的初始化频率拉开距离。弱通道失效是 MVMD 特有的现象。假如某个通道里有一个分量幅值特别小但它确实在物理上存在MVMD 的频率更新由于融入了所有通道的能量弱通道本身对频率估计的贡献不大。这是设计使然不算 bug但你会发现在结果里弱通道的这个模态幅值会被压缩。遇到这种情况我会在前处理阶段对通道做一定程度的幅度归一化让每个通道在频率估计里都有类似的发言权。边界效应则是所有频域分解方法都躲不开的问题。FFT 假设信号是周期重复的如果信号两端不在同一水平边界处的重构会出现比较明显的振荡。MVMD 对模态窄带性有约束边界振荡会被“摊”到附近若干采样点里。做法是分段加窗或者直接截掉两端各几个采样点做分析数据量充足时这是最省事的办法。5.3 通道预处理的优先级排序我给通道预处理排了一个优先级顺序。第一个是去直流。每个通道要减去均值否则第一模态容易吃进一个巨大的缓慢分量。第二个是坏通道剔除。如果某个通道完全被噪声覆盖或者出现持续性平直信号要果断在分解前删掉因为 MVMD 的频率更新会把这个通道的无效能量考虑进去把其他通道的频率估计带偏。第三个是通道归一化。各通道量纲差别大时比如一个通道是加速度单位、另一个是应变单位强制共享频率会造成弱通道被强通道压制。办法是让每个通道先除以其标准差做分解出结果后再乘回去。真实使用中预处理和数据本身质量的重要程度不亚于算法参数。MVMD 不是魔法输入脏数据输出的模态也只能是小了一点脏数据的“干净版本”。6. 这套算法适合哪些数据哪些别硬套6.1 适合的多变量信号场景我在实际里觉得 MVMD 最适合的场景是同步采集的传感器阵列数据、人体多导联生理信号、结构多测点振动响应以及旋转机械多个测点的故障振动信号。这些数据的共同特点是物理激励源是一致的所以各通道理应共享同一组模态频率只是幅值和相位随着测点位置变化。尤其是结构多测点的模态识别拿 MVMD 分解以后可以直接把每个测点在同一模态下的幅值和相位提出来做振型辨识。这和逐通道 VMD 之后再拼振型相比少了一步对齐少了一个误差源。6.2 不适合的场景与变通思路有两个场景不适合直接上 MVMD。第一个是各通道频率成分有明显频移的数据。比如非接触式测量里受多普勒效应影响每个通道同一物理模态的频率会不同。这种情况下硬要 MVMD 共享一个全局频率结果会是一个折中值各通道都拟合不好。第二个是通道数特别多、数据特别长的场景比如几千通道的阵列数据。MVMD 每轮迭代都要对 (K \times C) 个模态做频域更新计算量成倍上涨需要仔细做向量化否则内存和耗时都扛不住。遇到第二种场景可以先把通道分组每组里先跑 MVMD再用组级频率做后期融合。这是我对某些大规模数据常用的一种简化处理牺牲一点严格性换来工程可行性。MVMD 并不神秘它就是把 VMD 里的单通道信号换成了向量值信号并且让所有通道共享一个频率参数。理解这一点以后代码实现的核心就变得很直白模态更新的时候按通道独立做维纳滤波频率更新的时候把通道能量合起来算重心乘子更新则保持原来的 ADMM 节奏。我在实际调试中有一个体会值得分享出来初版代码不要追求极致的向量化先用清晰的三层循环写出来。三层循环慢但每个中间矩阵都能单独打印检查定位问题最快。模型调稳以后再做向量化把通道循环改成矩阵操作耗时能降下来但那是性能优化阶段该做的事不是理解算法阶段该做的事。另外如果某一步的结果和你预期完全对不上不要去怀疑数学推导先检查频谱索引和共轭对称那几行。九成的实现 bug 都藏在数组形状和频率轴的边界里。