基于希尔伯特-黄变换的信号处理方法研究

根据上面的定义，在理论上， h1(t) 满足：(1)极值点（极大值或极小值）数目与跨零点数目
相等或最多相差一个，(2)由局部极大值构成的上包络和由局部极小值构成的下包络的平均
值为零；即 h1(t) 应该是 IMF。实际上，由干包络线样条逼近的过冲和俯冲作用，会产生新
的极值影响原来极值的位置与大小；因此，分解得到的 h1(t) 并不完全满足 IMF 条件。用 h1(t)
(Tm ,Tn ) 和函数值 (U ,V ) 记为： Im = [Im (1), I m (2),L, Im (M )], I n = [I n (1), I n (2),L , I n (M )],
(3.14) (3.15)
Tm (i) = t Im ,U (i) = xIm , i = 1,L, M ,
针对非平稳数据的处理，许多学者通过修改傅立叶分析的全局表述，提出了一些改进措 施，如功率谱、小波分析、Wigner-Ville 分布、进化谱、经验正交函数分解(EOF)，这些都 基于傅立叶分析，所以它们仅限于线性系统。
N.E.Huang 等人于 1998 年提出了一种基于经验模式分解（EMD）的新数据分析方法—— 希尔伯特-黄变换，产生一个内在模式函数（IMF）集。这种分解是基于对与不同时间尺度 有关的能量的直接抽取。在 IMF 集的表达式中，它们都易于进行希尔伯特变换，而计算出 瞬时频率。因此我们可以关注时间轴和频率轴上的任何位置。从 IMF 的角度来说，这种分 解还可以看作是数据的扩展。然后，这些基于或源于数据的内在模式函数集可以做这种扩展 的基础，这种扩展由数据决定其是线性或者非线性，并且是完备的和几乎正交的。最重要的 是，它是自适应的。如同下文将要详细讨论的那样，局部性和自适应性是扩展非线性和非平 稳时间序列的必要条件；正交性对于我们对非平稳系统的基本选择来说，却不是一个必须的 标准。这种扩展的建设性理论是基于描述这种现象的振荡的实际时间尺度。通过希尔伯特变 换，源于 IMF 的局部能量和瞬时频率可以给我们一个完整的数据的能量－频率－时间分布。 这种表示称作希尔伯特谱，它对非线性和非平稳数据分析较为理想。
2. 希尔伯特-黄变换的理论
2.1 瞬时频率
我们已经接受了信号的瞬时能量或瞬时包络的概念，但对于瞬时频率的概念，却存在很 大的疑问。
在传统的傅立叶分析中，频率定义为在整个数据长度上有连续振幅的正弦或余弦函数。 作为这个定义的扩展，瞬时频率也需要与正弦或余弦函数关联。因此，我们需要至少一个正
-1-
http://www.paper.edu.cn
(2.1)
这个定义下的 X (t) 和Y (t ) 为复共轭对，所以我们可以得到解析信号 Z (t) ：
Z (t) = X (t) + iY (t) = a(t)eiθ (t) ,
(2.2)
这里
a(t) = [ X 2 (t) + Y 2 (t)]1 2 ,θ (t) = arctan( Y (t) ).
设时间序列信号为 X (t) ，它的上、下包络线分别为 u(t) 和 v(t) ，则上、下包络的平均
曲线 m(t) 为:
m(t) = 1 [u(t) + v(t)]
2
(3.1)
通过移动过程，用 X (t) 减去 m(t) 后剩余部分 h1(t) ，即:
h1(t) = X (t) − m(t),
(3.2)
信号的极大值和极小值数据集的规律, 预测附加的极值点，即所谓的端点延拓。
Baidu Nhomakorabea
设离散信号：
t ∈ [t(1),t(2),L,t(n)] = [t1,t2 ,L,tn ],
(3.12)
X (t) ∈ [x(t1), x(t2 ),L, x(tn )] = [x1, x2 ,L, xn ],
(3.13)
其采样步长为 ∆t , X (t) 有 M 个极大值和 N 个极小值，对应的序列下标 (Im , In ) 、时间
个 IMF, C1(t) 和信号的剩余部分为 r1(t) ，即:
C1(t) = hk (t),
(3.7)
r1(t) = X (t) − C1(t).
(3.8)
对信号的剩余部分 r1(t) 继续进行 EMD 分解，直到所得的剩余部分为一单调信号或其
值小于预先给定的值时，分解完毕。最终分解得到所有的 IMF 及余量：
symmetry of the upper and lower envelopes with respect to zero.
为了使用瞬时频率的这个特殊的定义，我们必须把一个随意数据集削减为 IMF 分量， 这样一个瞬时频率值可以被赋予各个 IMF 分量。从而，对于复杂数据，在一个局部时间我 们可以有多个瞬时频率。我们将引入经验模式分解方法来削减数据为所需的 IMF。
代替 X (t) ，与 h1(t) 相应的上、下包络线为 u1(t) 和 v1(t) ，重复移动过程，即:
m1(t) = [u1(t) + v1(t)] 2,
(3.3)
h2 (t) = h1(t) − m1(t),
(3.4)
…
[ ] mk −1(t) = uk −1(t) + vk −1(t) 2 ,
虑到样条插值的要求，所以只要在信号左、右两端分别延拓两个极大值和两个极小值即可[7]。
3个变量 a,ω,t ,互相相关，如3维平面的3条坐标轴。易得分别以 (ω, a), (ω, t), (t, a) 为底
的3维图像。其中任意一个3维图像都表示希尔伯特谱。而作为 (ω, t) 函数的幅度 a ，称为希
尔伯特幅度谱[4]， H (ω, t) 。
2.3.2 边际谱
同样，我们定义希尔伯特边际谱 h(ω) 如下[4]：
-2-
http://www.paper.edu.cn
2.3 希尔伯特谱和边际谱
2.3.1 希尔伯特谱
经过经验模式分解，我们得到信号的各个内在模式函数 C 1 , C2 ,...Cn , 有
X (t) = C1(t) + C2 (t) + LCn (t) + Rn (t)
(2.14)
根据式(2.1)- (2.4)，可得出所有的IMF的瞬时频率和瞬时幅度——均为时间 t 的函数。
弦或余弦波的全振动去定义局部的频率值。根据这个逻辑，任何短于一个整波的信号都没有 意义。对于非平稳数据而言，这种定义没有意义，因为其频率值时刻在变化。
对于任意一个时间序列 X (t) ，我们总是有它的 Hilbert[2]变换Y (t ) ：
∫ Y (t) = 1
∞ X (t ') dt ',
π −∞ t − t '
r2 = r1 − C2 ,L, rn = rn−1 − Cn ,
而原始的信号 X (t) 可表示为所有的 IMF 及余量之和[5]：
(3.9)
X (t) = C1(t) + C2 (t) + LCn (t) + Rn (t)
(3.10)
对于 IMF 的完备性和正交性，N.E.Huang 对此予以了详细的证明，本文不做过多描述[3]。
3
2
1
0
-1
-2
-3
80
100
120
140
160
180
200
220
Tim e(s )
图 2.1 典型的内在模式函数：相同数目的过零点和极值点，上下包络均值为 0。 Figure 2.1 A typical intrinsic mode function with the same number of zero crossings and extrema,and
(2.3)
通过 Hilbert 变换，瞬时频率定义为[3]:
X (t)
ω = dθ (t) .
(2.4)
dt
2.2 内在模式函数
实际上，定义一个有意义的瞬时频率的条件是函数相对于局部零均值是对称的，并且有 相等的过零点数和极值数。基于此，我们提出一种函数——内在模式函数（IMF）。
内在模式函数（IMF）满足两个条件：(1)在整个数据集中，极值数和过零数相等或最多 差 1；(2)在任意点，由局部极大值构成的包络和由局部极小值构成的包络的均值为零。图 2.1 给出了一个典型的 IMF。
(3.5)
hk (t) = hk −1 (t) − mk −1 (t),
(3.6)
-3-
http://www.paper.edu.cn
直到所得的 hk (t) 满足 IMF 条件:(1)极值点数目与过零点数目相等或最多相差一个，(2)由局
部极大值构成的上包络和由局部极小值构成的下包络的平均值趋近于零。这样就分解得第一
3.2 端点飞翼
3.2.1 问题的提出
对于原始信号进行 EMD 运算，首先就要得到其上、下包络，即根据信号极大极小值， 用三次样条曲线进行逼近。
由于所分析信号的有限长度、信号的两端点不能确定是极值，那么，在进行三次样条插 值的时候，必然使得信号的上下包络在信号的两端附近严重扭曲。这就是所谓的端点飞翼的 问题[6]。当处理低频或者多分量复杂信号时，会严重淹没信号的端部特征。
(3.16)
Tn (i) = tIn ,V (i) = xIn ,i = 1,L, N ,
(3.17)
在分解 IMF 时，端点延拓的目的是确保上、下包络都与端点相交，以便有与每一个信
号点相对应的局部平均值.而上、下包络是由极大值和极小值连结而成的，因此只要对极大
值和极小值进行延拓，而不必对信号本身进行延拓。极大值和极小值是相间分布的，同时考
从图 3.1 可以看出，对于多分量复杂信号的 EMD 分解，存在严重的端点飞翼效应。而 实际信号几乎都是复杂的多分量信号, 因此, 不能按照 EMD 的原始定义作 EMD 分解, 需要 研究合适的方法来抑制样条插值的端点飞翼, 同时又不能扭曲原始信号的端部特征。
3.2.2 解决方法：端点延拓
根据三次样条插值的特点, 必须在极大值和极小值数据集两端增加极大值和极小值点。 但是, 由于原始信号的两端点可能不是极值点, 必须进行合理的预测。为此, 基本的指导思 想就是在每次平滑过程中必须正确地确定添加极值点的位置和幅值。确定的方法是根据原始
http://www.paper.edu.cn
基于希尔伯特-黄变换的信号处理方法研究
龚程
北京邮电大学电信工程学院，北京（100876）
E-mail：handsome_gc@163.com
摘 要：针对目前非平稳信号处理方法效果不理想的情况，N.E.Huang 等人于 1998 年提出 了希尔伯特-黄变换(HHT)，其关键是利用经验模式分解(EMD)，将数据分解为确定的并且通 常幅度较小的内在模式函数（IMF），便于希尔伯特变换。本文首先介绍了希尔伯特-黄变 换的原理及铝电解的一些知识，然后将经验模式分解（EMD）应用于铝电解阳极效应检测。 关键词：非平稳时间序列，内在模式函数，经验模式分解，阳极效应 中图分类号：TP391
1. 问题的提出
数据分析是理论研究和实际应用中必不可少的一部分。由于数据的非理想化，因此，数 据分析通常有两个目的：一是确定构建模型所需的必要参数，二是确保我们所构建的模型尽 可能真实。然而，无论从物理测量还是从数学建模中得来的数据，都很可能存在如下一个或 多个问题：(a)时间跨度太短；(b)非平稳；(c)非线性。
∫ h(ω) = ∞ H (ω,t)dt −∞
(2.15)
其中 h(ω) 表明单位频率内的幅度分布（或者能量分布），代表着整个数据段幅度概率分布
的累加。
3. 希尔伯特-黄变换的实现
3.1 经验模式分解
与其他的信号处理方法相比，HHT 的创新点是引入了基于信号局部时间特性的 IMF， 以获得具有实际物理意义的瞬时频率。它主要由两部分组成，即 EMD 和 Hilbert 变换。EMD 通过多次的移动过程，一方面消除信号上的骑波（riding waves），另一方面对高低不平的 振幅进行平滑处理，使得每一个 IMF 具有如下两个特性：(1)极值点（极大值或极小值）数 目与过零点数目相等或最多相差一个，(2)由局部极大值构成的上包络和由局部极小值构成 的下包络的平均值为零。IMF 的上述两个特征，也即 EMD 分解结束的收敛准则。
历史上，傅立叶谱分析提供了一个检验全局能量－频率分布的一般方法，“谱”的概念几 乎等同于数据的傅里叶变换。自从引入傅里叶变换以来，它在数据分析上取得了巨大的成就， 而且被应用于各行各业。但是傅里叶变换有一些很严格的限制：(a)系统必须线性；(b)数据 必须严格周期平稳或宽平稳；否则，谱的结果没有任何实际意义[1]。