变分模态分解 (VMD) 详解
1. 引言与背景
变分模态分解(Variational Mode Decomposition,简称VMD)是由Konstantin Dragomiretskiy和Dominique Zosso在2014年提出的一种信号处理技术,旨在解决经验模态分解(Empirical Mode Decomposition,EMD)的一些固有缺点。信号分解是信号处理领域中的一个核心问题,其目标是将复杂的多组分信号分解为更简单、更有意义的组成部分(模态),以便于进一步分析和理解。在VMD出现之前,EMD是进行这类分解的主流方法。EMD由Huang等人在1990年代末提出,是一种算法性方法,通过递归地检测信号中的局部极值点,估计上下包络线,移除包络线平均值作为"低通"中心线,从而分离出高频振荡作为信号的"模态",然后在剩余的"低通"中心线上继续递归操作。虽然EMD在某些情况下确实能够将信号分解为主要模态,但由于其算法性质,EMD缺乏严格的数学理论支持,并且对噪声和采样敏感度较高。
近年来,研究人员提出了一些改进版本以解决EMD的这些问题。例如,更鲁棒的版本通过约束优化技术替代局部极值点的提取和插值,以形成包络线。其他方法采用部分变分方法,将信号显式建模为本征模态函数(IMF),或使用小波技术,如同步压缩(synchrosqueezing)和经验小波变换(Empirical Wavelet Transform,EWT)。然而,这些方法仍然存在一些局限性:大多数方法保留了EMD的递归筛选结构,无法进行后向误差修正;对噪声的处理能力有限;小波方法可能存在硬带限制;EWT需要预定义滤波器组边界。针对这些问题,VMD提出了一种全新的变分方法,能够自适应地确定相关频带,并同时估计相应的模态,从而在模态之间适当平衡误差。
2. 模态的定义与特性
在深入了解VMD之前,有必要明确模态的定义。在原始EMD中,模态被定义为局部极值点数量和零交叉点数量最多相差一个的信号。在后续工作中,这一定义稍微修改为所谓的本征模态函数(IMF),基于调制准则:
定义1(本征模态函数):本征模态函数是幅度调制-频率调制(AM-FM)信号,写作:
u ( t ) = A ( t ) cos ( ϕ ( t ) ) u(t) = A(t)\cos(\phi(t)) u(t)=A(t)cos(ϕ(t))
其中相位 ϕ ( t ) \phi(t) ϕ(t) 是非递减函数( ϕ ′ ( t ) > 0 \phi'(t) > 0 ϕ′(t)>0),包络 A ( t ) A(t) A(t) 非负( A ( t ) ≥ 0 A(t) \ge 0 A(t)≥0),并且非常重要的是,包络和瞬时频率 ω ( t ) = ϕ ′ ( t ) \omega(t) = \phi'(t) ω(t)=ϕ′(t) 的变化比相位变化慢得多,即 A ′ ( t ) ≪ A ( t ) ϕ ′ ( t ) A'(t) \ll A(t)\phi'(t) A′(t)≪A(t)ϕ′(t)和 ϕ ′ ′ ( t ) ≪ ( ϕ ′ ( t ) ) 2 \phi''(t) \ll (\phi'(t))^2 ϕ′′(t)≪(ϕ′(t))2。
这意味着在足够长的时间间隔 [ t − ϵ , t + ϵ ] [t-\epsilon, t+\epsilon] [t−ϵ,t+ϵ] 内,模态可以被视为具有振幅 A ( t ) A(t) A(t) 和瞬时频率 ω ( t ) = ϕ ′ ( t ) \omega(t) = \phi'(t) ω(t)=ϕ′(t) 的纯谐波信号。需要注意的是,较新的IMF定义比原始定义更为严格:虽然所有符合上述IMF定义的信号都满足原始EMD模态的性质,但反之并不一定成立。
新的IMF定义的直接结果是带宽受限,这是VMD模型中进行模态分离的核心假设。如果 ω k \omega_k ωk 是模态的平均频率,则其实际带宽受到两个因素的影响:
- 瞬时频率相对于其中心(载波频率 f c f_c fc)的最大偏差 Δ f \Delta f Δf
- 这种偏差的变化率 f F M f_{FM} fFM
根据卡森规则(Carson’s rule),频率调制信号的带宽可以估计为:
B W F M = 2 ( Δ f + f F M ) BW_{FM} = 2(\Delta f + f_{FM}) BWFM=2(Δf+fFM)
此外,调制FM信号幅度的包络本身的带宽(由其最高频率 f A M f_{AM} fAM 给出)会进一步拓宽频谱:
定义2(总实用IMF带宽):IMF的总实用带宽估计为:
B W A M − F M = 2 ( Δ f + f F M + f A M ) BW_{AM-FM} = 2(\Delta f + f_{FM} + f_{AM}) BWAM−FM=2(Δf+fFM+fAM)
根据具体的IMF,这些项中的任何一个都可能占主导地位。
3. 信号处理工具基础
为了理解VMD的工作原理,需要回顾几个信号处理中的基本概念,这些概念构成了VMD变分模型的构建基块。
3.1 维纳滤波
首先考虑一个简单的去噪问题。假设观测信号 f f f 是原始信号受到加性零均值高斯噪声影响的版本:
f ( t ) = s ( t ) + n ( t ) f(t) = s(t) + n(t) f(t)=s(t)+n(t)
恢复未知信号 s s s 是一个典型的不适定逆问题,通常使用Tikhonov正则化处理:
min s ∥ f − s ∥ 2 2 + α ∥ L s ∥ 2 2 \min_s \|f - s\|_2^2 + \alpha \|Ls\|_2^2 smin∥f−s∥22+α∥Ls∥22
其中 L L L 是某种微分算子。这个问题的欧拉-拉格朗日方程在傅里叶域中很容易求解:
s ^ ( ω ) = f ^ ( ω ) 1 + α ∣ L ^ ( ω ) ∣ 2 \hat{s}(\omega) = \frac{\hat{f}(\omega)}{1 + \alpha|\hat{L}(\omega)|^2} s^(ω)=1+α∣L^(ω)∣2f^(ω)
其中 s ^ \hat{s} s^、 f ^ \hat{f} f^ 和 L ^ \hat{L} L^ 分别是 s s s、 f f f 和 L L L 的傅里叶变换。显然,恢复的信号是输入信号围绕 ω = 0 \omega = 0 ω=0 的低通窄带选择。这个解对应于与维纳滤波器的卷积,其中 α \alpha α 表示白噪声的方差,而信号具有低通功率谱先验。
3.2 希尔伯特变换与解析信号
希尔伯特变换是一种线性、平移不变的算子,将所有一维余弦函数映射到相应的正弦函数。它是一个全通滤波器,其传递函数为 H ^ ( ω ) = − j ⋅ sgn ( ω ) \hat{H}(\omega) = -j \cdot \text{sgn}(\omega) H^(ω)=−j⋅sgn(ω)。希尔伯特变换的最重要用途是从纯实信号构造解析信号,这是由Gabor提出的:
定义4(解析信号):设 u ( t ) u(t) u(t) 是一个纯实信号。复值解析信号定义为:
A ( t ) = u ( t ) + j ⋅ H { u } ( t ) A(t) = u(t) + j \cdot \mathcal{H}\{u\}(t) A(t)=u(t)+j⋅H{u}(t)
这个解析信号有几个重要特性:
-
复指数项 e j ϕ ( t ) e^{j\phi(t)} ejϕ(t) 是一个相量,描述复信号随时间的旋转, ϕ ( t ) \phi(t) ϕ(t) 是相位,而振幅由实包络 A ( t ) A(t) A(t) 控制。这种表示在分析时变振幅和瞬时频率时特别有用,后者定义为 ω ( t ) = ϕ ′ ( t ) \omega(t) = \phi'(t) ω(t)=ϕ′(t)。
-
解析信号的频谱是单边的,只包含非负频率。
-
从解析信号中,原始(实)信号可以简单地作为实部恢复: u ( t ) = Re { A ( t ) } u(t) = \text{Re}\{A(t)\} u(t)=Re{A(t)}。
3.3 频率混合与外差解调
最后一个需要回顾的概念是频率混合原理。混合是将两个信号非线性组合的过程,从而在输出中引入交叉频率项。最简单的混频器是乘法。将两个频率分别为 ω 1 \omega_1 ω1 和 ω 2 \omega_2 ω2 的实信号相乘,会在输出中产生频率为 ω 1 + ω 2 \omega_1 + \omega_2 ω1+ω2 和 ω 1 − ω 2 \omega_1 - \omega_2 ω1−ω2 的混合频率。
在VMD中,作者混合了两个相应的解析信号:
A 1 ( t ) ⋅ e − j ω 2 t = u 1 ( t ) e j ϕ 1 ( t ) e − j ω 2 t = u 1 ( t ) e j ( ϕ 1 ( t ) − ω 2 t ) A_1(t) \cdot e^{-j\omega_2 t} = u_1(t)e^{j\phi_1(t)}e^{-j\omega_2 t} = u_1(t)e^{j(\phi_1(t)-\omega_2 t)} A1(t)⋅e−jω2t=u1(t)ejϕ1(t)e−jω2t=u1(t)ej(ϕ1(t)−ω2t)
因此,混合信号自动变为"单音"(仅由单一频率构成)。在傅里叶术语中,这对应于以下变换对:
e j ω 0 t ↔ 2 π δ ( ω − ω 0 ) e^{j\omega_0 t} \leftrightarrow 2\pi\delta(\omega-\omega_0) ejω0t↔2πδ(ω−ω0)
因此,将解析信号与纯指数相乘会导致简单的频率移位。
4. 变分模态分解模型
VMD的目标是将实值输入信号 f f f 分解为离散数量的子信号(模态) { u k } \{u_k\} {uk},这些子信号具有特定的稀疏性质,并且能够重构输入信号。每个模态的稀疏性先验被选为其在频谱域中的带宽。换句话说,假设每个模态主要集中在一个待定的中心脉动 ω k \omega_k ωk 周围。
为了评估模态的带宽,VMD采用以下方案:
- 对于每个模态 u k u_k uk,通过希尔伯特变换计算相关的解析信号,以获得单边频谱。
- 对于每个模态,通过与调谐到各自估计中心频率的复指数混合,将模态的频谱移到"基带"。
- 通过解调信号的高斯平滑度(即梯度的平方 L 2 L^2 L2 范数)估计带宽。
由此得到的约束变分问题表述为:
min { u k } , { ω k } ∑ k ∥ ∂ t [ ( δ ( t ) + j π t ) ∗ u k ( t ) ] e − j ω k t ∥ 2 2 \min_{\{u_k\},\{\omega_k\}} \sum_{k} \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 {uk},{ωk}mink∑ ∂t[(δ(t)+πtj)∗uk(t)]e−jωkt 22
s.t. ∑ k u k = f \text{s.t.} \sum_{k} u_k = f s.t.k∑uk=f
其中 { u k } \{u_k\} {uk} 和 { ω k } \{\omega_k\} {ωk} 分别是所有模态及其中心频率的集合的简写, ∗ * ∗ 表示卷积, δ ( t ) \delta(t) δ(t) 是狄拉克脉冲函数。
为了解决重建约束,作者使用二次惩罚项和拉格朗日乘子 λ \lambda λ,使问题成为无约束问题。因此,增广拉格朗日函数为:
L ( { u k } , { ω k } , λ ) = α ∑ k ∥ ∂ t [ ( δ ( t ) + j π t ) ∗ u k ( t ) ] e − j ω k t ∥ 2 2 + ∥ f ( t ) − ∑ k u k ( t ) ∥ 2 2 + ⟨ λ ( t ) , f ( t ) − ∑ k u k ( t ) ⟩ \mathcal{L}(\{u_k\}, \{\omega_k\}, \lambda) = \alpha \sum_{k} \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 + \left\| f(t) - \sum_{k} u_k(t) \right\|_2^2 + \left\langle \lambda(t), f(t) - \sum_{k} u_k(t) \right\rangle L({uk},{ωk},λ)=αk∑ ∂t[(δ(t)+πtj)∗uk(t)]e−jωkt 22+ f(t)−k∑uk(t) 22+⟨λ(t),f(t)−k∑uk(t)⟩
原始最小化问题的解决方案现在是增广拉格朗日的鞍点,在一系列称为乘子交替方向法(ADMM)的迭代子优化中找到。
4.1 ADMM优化算法
ADMM是一种解决带约束优化问题的有效方法,特别是当目标函数可以分解为由不同变量优化的子问题时。在VMD中,优化过程如下:
算法1:VMD的ADMM优化概念
初始化
{
u
k
1
}
\{u_k^1\}
{uk1},
{
ω
k
1
}
\{\omega_k^1\}
{ωk1},
λ
1
\lambda^1
λ1,
n
←
0
n\leftarrow 0
n←0
重复
n
←
n
+
1
n\leftarrow n+1
n←n+1
对于每个
k
k
k 做
更新
u
k
u_k
uk:
u
k
n
+
1
←
arg
min
u
k
L
(
{
u
i
n
+
1
}
i
<
k
,
u
k
,
{
u
i
n
}
i
>
k
,
{
ω
i
n
}
,
λ
n
)
u_k^{n+1} \leftarrow \arg\min_{u_k} \mathcal{L}(\{u_i^{n+1}\}_{i<k}, u_k, \{u_i^n\}_{i>k}, \{\omega_i^n\}, \lambda^n)
ukn+1←argminukL({uin+1}i<k,uk,{uin}i>k,{ωin},λn)
结束
对于每个
k
k
k 做
更新
ω
k
\omega_k
ωk:
ω
k
n
+
1
←
arg
min
ω
k
L
(
{
u
i
n
+
1
}
,
{
ω
i
n
+
1
}
i
<
k
,
ω
k
,
{
ω
i
n
}
i
>
k
,
λ
n
)
\omega_k^{n+1} \leftarrow \arg\min_{\omega_k} \mathcal{L}(\{u_i^{n+1}\}, \{\omega_i^{n+1}\}_{i<k}, \omega_k, \{\omega_i^n\}_{i>k}, \lambda^n)
ωkn+1←argminωkL({uin+1},{ωin+1}i<k,ωk,{ωin}i>k,λn)
结束
双重上升:
λ
n
+
1
←
λ
n
+
τ
(
f
−
∑
k
u
k
n
+
1
)
\lambda^{n+1} \leftarrow \lambda^n + \tau(f - \sum_k u_k^{n+1})
λn+1←λn+τ(f−∑kukn+1)
直到收敛:
∑
k
∥
u
k
n
+
1
−
u
k
n
∥
2
2
/
∥
u
k
n
∥
2
2
<
ϵ
\sum_k \|u_k^{n+1} - u_k^n\|_2^2 / \|u_k^n\|_2^2 < \epsilon
∑k∥ukn+1−ukn∥22/∥ukn∥22<ϵ
4.2 关于 u k u_k uk 的最小化
对于更新模态 u k u_k uk,子问题可以重写为:
min u k α ∥ ∂ t [ ( δ ( t ) + j π t ) ∗ u k ( t ) ] e − j ω k t ∥ 2 2 + ∥ f ( t ) − ∑ i u i ( t ) + λ ( t ) 2 ∥ 2 2 \min_{u_k} \alpha \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 + \left\| f(t) - \sum_{i} u_i(t) + \frac{\lambda(t)}{2} \right\|_2^2 ukminα ∂t[(δ(t)+πtj)∗uk(t)]e−jωkt 22+ f(t)−i∑ui(t)+2λ(t) 22
利用 L 2 L^2 L2 范数下的 Parseval/Plancherel 傅里叶等距性,这个问题可以在频谱域中解决:
min u k α ∫ ∣ ω A ^ k ( ω − ω k ) ∣ 2 d ω + ∫ ∣ f ^ ( ω ) − ∑ i u ^ i ( ω ) + λ ^ ( ω ) 2 ∣ 2 d ω \min_{u_k} \alpha \int |\omega \hat{A}_k(\omega-\omega_k)|^2 d\omega + \int \left|\hat{f}(\omega) - \sum_i \hat{u}_i(\omega) + \frac{\hat{\lambda}(\omega)}{2}\right|^2 d\omega ukminα∫∣ωA^k(ω−ωk)∣2dω+∫ f^(ω)−i∑u^i(ω)+2λ^(ω) 2dω
通过变量变换和利用实信号在重建保真度项中的厄米对称性,最优解为:
u ^ k ( ω ) = f ^ ( ω ) − ∑ i ≠ k u ^ i ( ω ) + λ ^ ( ω ) 2 1 + 2 α ( ω − ω k ) 2 \hat{u}_k(\omega) = \frac{\hat{f}(\omega) - \sum_{i \neq k} \hat{u}_i(\omega) + \frac{\hat{\lambda}(\omega)}{2}}{1 + 2\alpha(\omega - \omega_k)^2} u^k(ω)=1+2α(ω−ωk)2f^(ω)−∑i=ku^i(ω)+2λ^(ω)
这是当前残差的维纳滤波,信号先验为 1 ( ω − ω k ) 2 \frac{1}{(\omega - \omega_k)^2} (ω−ωk)21。模态在时间域中则是这种滤波解析信号的实部。
4.3 关于 ω k \omega_k ωk 的最小化
中心频率的更新只涉及带宽先验项:
min ω k ∥ ∂ t [ ( δ ( t ) + j π t ) ∗ u k ( t ) ] e − j ω k t ∥ 2 2 \min_{\omega_k} \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 ωkmin ∂t[(δ(t)+πtj)∗uk(t)]e−jωkt 22
这个问题在傅里叶域中优化,解为:
ω k = ∫ 0 ∞ ω ∣ u ^ k ( ω ) ∣ 2 d ω ∫ 0 ∞ ∣ u ^ k ( ω ) ∣ 2 d ω \omega_k = \frac{\int_0^{\infty} \omega |\hat{u}_k(\omega)|^2 d\omega}{\int_0^{\infty} |\hat{u}_k(\omega)|^2 d\omega} ωk=∫0∞∣u^k(ω)∣2dω∫0∞ω∣u^k(ω)∣2dω
这将新的 ω k \omega_k ωk 设置为对应模态功率谱的重心。
4.4 完整的VMD算法
将所有子优化的解决方案结合起来,并直接在适当的地方在傅里叶域中优化,得到完整的变分模态分解算法:
算法2:VMD的完整优化
初始化
{
u
^
k
1
}
\{\hat{u}_k^1\}
{u^k1},
{
ω
k
1
}
\{\omega_k^1\}
{ωk1},
λ
^
1
\hat{\lambda}^1
λ^1,
n
←
0
n\leftarrow 0
n←0
重复
n
←
n
+
1
n\leftarrow n+1
n←n+1
对于每个
k
k
k 做
更新所有
ω
\omega
ω 的
u
^
k
\hat{u}_k
u^k:
u
^
k
n
+
1
(
ω
)
←
f
^
(
ω
)
−
∑
i
≠
k
u
^
i
n
(
ω
)
+
λ
^
n
(
ω
)
2
1
+
2
α
(
ω
−
ω
k
n
)
2
\hat{u}_k^{n+1}(\omega) \leftarrow \frac{\hat{f}(\omega) - \sum_{i \neq k} \hat{u}_i^n(\omega) + \frac{\hat{\lambda}^n(\omega)}{2}}{1 + 2\alpha(\omega - \omega_k^n)^2}
u^kn+1(ω)←1+2α(ω−ωkn)2f^(ω)−∑i=ku^in(ω)+2λ^n(ω)
更新
ω
k
\omega_k
ωk:
ω
k
n
+
1
←
∫
0
∞
ω
∣
u
^
k
n
+
1
(
ω
)
∣
2
d
ω
∫
0
∞
∣
u
^
k
n
+
1
(
ω
)
∣
2
d
ω
\omega_k^{n+1} \leftarrow \frac{\int_0^{\infty} \omega |\hat{u}_k^{n+1}(\omega)|^2 d\omega}{\int_0^{\infty} |\hat{u}_k^{n+1}(\omega)|^2 d\omega}
ωkn+1←∫0∞∣u^kn+1(ω)∣2dω∫0∞ω∣u^kn+1(ω)∣2dω
结束
双重上升,对所有
ω
\omega
ω 做:
λ
^
n
+
1
(
ω
)
←
λ
^
n
(
ω
)
+
τ
(
f
^
(
ω
)
−
∑
k
u
^
k
n
+
1
(
ω
)
)
\hat{\lambda}^{n+1}(\omega) \leftarrow \hat{\lambda}^n(\omega) + \tau\left(\hat{f}(\omega) - \sum_k \hat{u}_k^{n+1}(\omega)\right)
λ^n+1(ω)←λ^n(ω)+τ(f^(ω)−∑ku^kn+1(ω))
直到收敛:
∑
k
∥
u
^
k
n
+
1
−
u
^
k
n
∥
2
2
/
∥
u
^
k
n
∥
2
2
<
ϵ
\sum_k \|\hat{u}_k^{n+1} - \hat{u}_k^n\|_2^2 / \|\hat{u}_k^n\|_2^2 < \epsilon
∑k∥u^kn+1−u^kn∥22/∥u^kn∥22<ϵ
4.5 精确重构与去噪
在VMD中,拉格朗日乘子的作用是强制执行约束,而二次惩罚则改善收敛性。如果精确重构不是目标,特别是在存在(强)噪声且不应将其包含在分解中的情况下,仅使用二次惩罚而省略拉格朗日乘子是合适的选择。实际上,单独的二次惩罚代表了与加性高斯噪声相关的最小二乘数据保真度先验。拉格朗日乘子可以通过选择其更新参数 τ = 0 \tau = 0 τ=0 而有效地保持为0,从而被关闭。
4.6 边界和周期性
在信号处理中,我们通常处理时间和分辨率都有限的信号。当考虑短时间信号时,隐含的假设是所考虑的信号只是无限长周期信号的一个周期提取。因此,短区间上的"一般趋势"函数的频谱包含大量高频谐波,因为我们实际上正在看周期性锯齿函数的频谱。在时间域中,周期化函数在域的端点处是不连续的,这严重影响了平滑项。为了解决这一问题,最简单的方法是对信号进行简单的镜像扩展,在每侧延长一半长度,这有效地对应于纽曼边界条件(即时间域边界处的导数消失)。
5. 实验结果与分析
作者进行了一系列实验来评估VMD的有效性,特别是与EMD相比较。
5.1 音调与采样
当输入信号是纯谐波时,模态分解应该准确地输出该谐波。然而,EMD在频率增加时表现出重要的抖动,并且当采样频率是音调频率的偶数倍时,存在近乎完美重构的尖峰。VMD的相对重构误差在很大程度上独立于谐波的频率,并且由收敛容差水平 ϵ \epsilon ϵ 很好地控制。只有在接近奈奎斯特频率的频率处,重构误差才会增加。
5.2 音调分离
当输入信号由两个纯谐波组成时,EMD在某些频率比例下难以正确分离它们,特别是当两个频率相近时。相比之下,VMD几乎在整个域上都能实现良好的音调分离,除了接近奈奎斯特频率的情况。特别是,对于频率接近的谐波,VMD的分解质量并没有显著恶化。
5.3 抗噪声性
为了测试VMD对噪声的鲁棒性,作者使用了一个受噪声影响的三谐波信号:
f n ( t ) = cos ( 4 π t ) + 1 4 cos ( 48 π t ) + 1 16 cos ( 576 π t ) + η f_n(t) = \cos(4\pi t) + \frac{1}{4}\cos(48\pi t) + \frac{1}{16}\cos(576\pi t) + \eta fn(t)=cos(4πt)+41cos(48πt)+161cos(576πt)+η
其中 η ∼ N ( 0 , σ ) \eta \sim \mathcal{N}(0,\sigma) η∼N(0,σ) 是高斯加性噪声。对于 σ = 0.1 \sigma = 0.1 σ=0.1(相对于最强和最弱谐波的振幅相当显著),VMD仍然能够相当准确地恢复最强的低频信号和中等强度的中频信号。对于弱的高频信号,VMD正确地调整了第三个中心频率,但恢复的模态受到噪声的严重影响。相比之下,EMD在相同信号上产生7个估计模态,但没有一个模态对应于纯谐波,并且所有模态都包含显著的噪声成分。
5.4 初始化与收敛
VMD算法不保证收敛到全局最小值,因此结果取决于初始化,特别是成功率取决于噪声水平。对于足够低的噪声水平,拉格朗日乘子可以强制精确重构输入信号,从而即使在随机初始化的情况下也能实现100%的成功率。随着噪声水平的增加,拉格朗日乘子变得不那么有用,因为此时希望从估计的模态中消除噪声。有趣的是,数据中的一些噪声对于找到正确的中心频率是有益的,类似于模拟退火中的非零温度,能够帮助算法逃离错误的局部最小值。
5.5 模态数量的选择
VMD需要预先指定模态的数量 K K K,这是许多聚类和分割算法的常见限制。当指定的模态数量太少时(欠分割),一个模态可能被相邻模态共享,或者主要被丢弃。当指定的模态数量太多时(过分割),额外的模态可能主要包含噪声,或者出现模态复制现象,其中多个模态具有相同的中心频率。信号分解为多个AM-FM信号并不是唯一的。例如,两个频率接近的谐波可以被感知为一个中间谐波的振幅调制,而不是两个单独的谐波。VMD可以根据指定的模态数量 K K K、噪声水平和带宽限制 α \alpha α 被引导为将这样的信号识别为两个单独的谐波或一个单一的AM-FM模态。
5.6 复杂信号的分解
作者还测试了VMD在多种复杂信号上的性能,包括:
-
具有线性趋势和谐波的信号:VMD能够将频谱很好地分割为不同的模态,每个模态在其各自的中心频率周围占主导地位。与EMD相比,VMD的分解更加干净,线性趋势和中频信号中的伪振荡更少。
-
具有二次趋势、啁啾信号和分段常数频率的信号:VMD估计的中心频率精确地收敛到预期频率,并且频谱分割清晰可见。虽然EMD在恢复二次趋势方面表现更好,但无法分离分段常数频率信号的两部分。
-
具有内波频率调制的信号:尽管这种信号违反了窄带假设,但VMD仍然能够捕获主要频率,虽然高阶谐波在两个模态之间共享。
-
锯齿信号:虽然锯齿信号的高阶谐波遍布整个频谱,违反了VMD模态的有限带宽假设,但VMD仍然能够捕获相关的中心频率,并实现相对良好的信号分离。
-
心电图(ECG)实际信号:VMD能够将ECG信号分解为多个模态,包括低频基线振荡、心跳频率的主要ECG模态,以及围绕非正弦尖峰的高阶波包。通过丢弃低频基线振荡和大部分高频噪声,可以重构出"干净"的ECG信号。
6. VMD的优势与局限性
6.1 优势
-
数学基础:相比EMD的算法性质,VMD具有坚实的变分理论基础。
-
并行提取:模态是并行估计的,而不是像EMD那样递归筛选,这允许适当平衡模态之间的误差。
-
抗噪声性:VMD对噪声的鲁棒性强,特别是当使用无拉格朗日乘子的去噪模式时。
-
频率分离:VMD能够有效分离频率接近的组件,这是EMD的一个主要限制。
-
自适应频带:VMD自动确定相关的频带,而不需要预定义滤波器组边界。
6.2 局限性
-
固定模态数量:VMD需要预先指定模态的数量,这在实际应用中可能需要先验知识。
-
边界效应:由于使用基于 L 2 L^2 L2 的平滑项,VMD对域边界和内部的跳跃过度惩罚,导致边界效应。
-
全局趋势处理:VMD不显式处理信号中的全局趋势,虽然可以通过固定一个模态的中心频率为0来部分解决这个问题。
-
非平稳信号:VMD不直接适用于非平稳信号,如长期EEG记录,因为模态的频谱带在时间上可能变化很大,并且在全局上重叠。
-
初始化敏感性:VMD的结果取决于初始化,特别是在噪声较大的情况下。
7. 总结与展望
VMD通过变分框架解决了信号分解问题,将信号表示为具有有限带宽的本征模态函数的集合。在这个框架中,每个模态的带宽被定义为希尔伯特补充解析信号(经过移频到基带)的 L 2 L^2 L2 范数。优化问题通过ADMM方法高效解决,其中模态通过维纳滤波更新,中心频率通过重心计算确定,并且拉格朗日乘子通过双重上升强制执行精确信号重构。与EMD和其他方法相比,VMD在音调检测、音调分离和抗噪声性方面表现出色。它能够处理各种复杂信号,从简单的谐波组合到实际的ECG数据,并提供有意义的分解。此外,VMD与维纳滤波有着密切的联系,这表明在处理噪声方面具有某种最优性。
尽管VMD具有许多优势,但它也有一些局限性。为了进一步发展VMD,可以考虑以下方向:
-
自动确定模态数量:开发能够自动确定最佳模态数量的方法,可能基于信息理论准则或贝叶斯模型选择。
-
改进边界处理:采用更先进的边界处理技术,减少边界效应对分解的影响。
-
扩展到非平稳信号:通过窗口化方法或自适应模型,使VMD能够处理频谱特性随时间变化的非平稳信号。
-
高维扩展:将VMD扩展到二维或更高维度的信号,如图像或体积数据。
-
集成学习:将VMD与机器学习方法结合,处理更复杂、更大规模的实际数据集。
8. VMD的数学原理深入解析
8.1 变分框架的理论基础
VMD的核心是将信号分解问题表述为变分问题。变分方法在数学物理中有着深厚的历史背景,特别是在最小作用原理和求解偏微分方程方面。在信号处理中,变分方法提供了一种强大的框架,可以将先验知识(如稀疏性、平滑性等)纳入信号恢复问题。VMD的变分框架建立在以下假设上:每个模态都具有有限的带宽,集中在其各自的中心频率周围。这种假设与本征模态函数(IMF)的定义一致,后者是幅度调制-频率调制(AM-FM)信号,其包络和瞬时频率变化比相位变化慢得多。
8.2 从希尔伯特变换到解析信号
希尔伯特变换在VMD中起着关键作用,用于构造解析信号。希尔伯特变换的数学定义是:
H { u } ( t ) = 1 π p.v. ∫ − ∞ ∞ u ( τ ) t − τ d τ \mathcal{H}\{u\}(t) = \frac{1}{\pi} \text{p.v.} \int_{-\infty}^{\infty} \frac{u(\tau)}{t - \tau} d\tau H{u}(t)=π1p.v.∫−∞∞t−τu(τ)dτ
其中p.v.表示柯西主值。在频域中,希尔伯特变换对应于传递函数 H ^ ( ω ) = − j ⋅ sgn ( ω ) \hat{H}(\omega) = -j \cdot \text{sgn}(\omega) H^(ω)=−j⋅sgn(ω),这意味着它将信号的每个频率分量的相位移动 − π / 2 -\pi/2 −π/2,同时保持其振幅不变。
解析信号 A ( t ) = u ( t ) + j ⋅ H { u } ( t ) A(t) = u(t) + j \cdot \mathcal{H}\{u\}(t) A(t)=u(t)+j⋅H{u}(t) 是一个复值信号,其频谱仅包含非负频率。这使得解析信号特别适合用于分析信号的瞬时振幅和频率特性。对于满足Bedrosian定理条件的IMF信号 u ( t ) = A ( t ) cos ( ϕ ( t ) ) u(t) = A(t)\cos(\phi(t)) u(t)=A(t)cos(ϕ(t)),其解析信号直接继承相同的振幅函数: A ( t ) = A ( t ) e j ϕ ( t ) A(t) = A(t)e^{j\phi(t)} A(t)=A(t)ejϕ(t)。
8.3 维纳滤波与VMD的关系
维纳滤波是一种最优线性估计器,用于从噪声信号中恢复期望信号。其公式为:
s ^ ( ω ) = P s ( ω ) P s ( ω ) + P n ( ω ) f ^ ( ω ) \hat{s}(\omega) = \frac{P_s(\omega)}{P_s(\omega) + P_n(\omega)} \hat{f}(\omega) s^(ω)=Ps(ω)+Pn(ω)Ps(ω)f^(ω)
其中 P s ( ω ) P_s(\omega) Ps(ω) 和 P n ( ω ) P_n(\omega) Pn(ω) 分别是信号和噪声的功率谱密度。
在VMD中,模态更新公式为:
u ^ k ( ω ) = f ^ ( ω ) − ∑ i ≠ k u ^ i ( ω ) + λ ^ ( ω ) 2 1 + 2 α ( ω − ω k ) 2 \hat{u}_k(\omega) = \frac{\hat{f}(\omega) - \sum_{i \neq k} \hat{u}_i(\omega) + \frac{\hat{\lambda}(\omega)}{2}}{1 + 2\alpha(\omega - \omega_k)^2} u^k(ω)=1+2α(ω−ωk)2f^(ω)−∑i=ku^i(ω)+2λ^(ω)
这可以解释为对残差信号进行维纳滤波,其中信号先验为 1 ( ω − ω k ) 2 \frac{1}{(\omega - \omega_k)^2} (ω−ωk)21,这对应于窄带信号围绕中心频率 ω k \omega_k ωk。这表明VMD可以看作是将经典维纳滤波器推广到多个自适应频带的过程。
8.4 ADMM优化的收敛性分析
乘子交替方向法(ADMM)是一种有效解决带约束优化问题的算法,特别是当目标函数可以分解为由不同变量优化的子问题时。ADMM的收敛性已在许多情况下得到证明,特别是当目标函数是凸的时。
对于VMD,虽然整体问题不是凸的(由于模态和中心频率的耦合非线性),但每个子问题(固定 ω k \omega_k ωk 时关于 u k u_k uk 的优化,以及固定 u k u_k uk 时关于 ω k \omega_k ωk 的优化)都是凸的。这种结构使得ADMM特别适合用于VMD问题。
尽管全局收敛性无法保证(由于非凸性),但实际经验表明,ADMM在大多数情况下都能收敛到实用的局部最优解。收敛速度取决于参数 α \alpha α(控制带宽约束的强度)和 τ \tau τ(控制拉格朗日乘子更新的步长)的选择,以及初始化条件。
9. VMD的实现与工程考虑
9.1 参数选择
VMD的性能在很大程度上取决于几个关键参数的选择:
-
模态数量K:这应该基于先验知识或通过实验确定。过少的模态会导致欠拟合,而过多的模态可能导致过拟合或模态复制。
-
带宽约束α:较大的α值强制模态的带宽更窄,而较小的α允许更宽的带宽。这个参数应该根据信号的预期特性来调整。
-
拉格朗日乘子更新步长τ:这影响拉格朗日乘子的更新速度和整体收敛性。太大的τ可能导致不稳定,而太小的τ可能导致收敛缓慢。
-
初始化:模态和中心频率的初始值会影响算法的收敛和最终结果。在实践中,经常使用随机初始化或基于信号频谱的初始猜测。
9.2 计算复杂度
VMD的主要计算成本来自于每次迭代中模态在频域的更新,这涉及傅里叶变换和逆变换。对于长度为N的信号和K个模态,每次迭代的复杂度为O(KN log N)。总的计算成本还取决于迭代次数,这又受到收敛条件和参数选择的影响。虽然VMD的计算复杂度高于某些简单的滤波方法,但它通常低于EMD,特别是对于复杂信号,并且可以通过并行计算和优化实现进一步降低。
9.3 实际实现中的考虑
在实际实现VMD时,需要注意以下几点:
-
边界处理:如前所述,镜像扩展是一种简单且有效的方法来减少边界效应。
-
数值稳定性:在计算频域更新和中心频率估计时,应注意避免除以零或接近零的值,特别是当信号功率很小时。
-
收敛准则:选择适当的收敛准则对于确保算法停止在合理的结果上很重要。常用的标准是模态之间的相对变化小于预定阈值。
-
并行计算:VMD算法中的模态更新是独立的,可以在多核处理器或GPU上并行执行,从而加速计算。
参考文献
- K. Dragomiretskiy and D. Zosso, “Variational Mode Decomposition,” in IEEE Transactions on Signal Processing, vol. 62, no. 3, pp. 531-544, Feb.1, 2014.