基于双重对称延拓的HHT端点效应抑制方法

项洪，吴琛, 杜喜朋

(1.福建工程学院 土木工程学院，福建 福州 350118; 2.福建省土木工程新技术与信息化重点实验室，福建 福州 350118)



基于双重对称延拓的HHT端点效应抑制方法

项洪1,2，吴琛1,2, 杜喜朋1,2

(1.福建工程学院 土木工程学院，福建 福州 350118; 2.福建省土木工程新技术与信息化重点实验室，福建 福州 350118)

摘要：经验模态分解( empirical mode decomposition，EMD)存在端点飞翼的固有问题，使信号两端出现扭曲失真。为抑制EMD端点效应，文章提出双重对称延拓法，以端部数据对称延拓作为首次延拓，以极值点对称延拓作为二次延拓，该方法可同时实现EMD分解和Hilbert变换两阶段端点效应的抑制。仿真信号和地震响应的工程算例通过IMF分量对比、瞬时频率对比，以及整体正交性、相关系数等指标的对比，验证了该方法具有良好的端点效应抑制效果,同时还能改善IMF分量之间的正交性，并节约计算时间。

关键词：Hilbert-Huang变换; 经验模态分解; 端点效应; 双重对称延拓; 评价指标

Hilbert-Huang变换[1](Hilbert-Huang transform , HHT) 包括经验模态分解(empirical mode decomposition, EMD)和Hilbert谱分析两个阶段。HHT具有完全自适应性的优点，处理非线性、非平稳的信号具有明显优势；自1998年Norden E.Huang等人提出后，被广泛应用于机械故障诊断[2-3]、地震工程[4-5]等领域。但在实际工程应用中，HHT还存在一些问题：如端点问题、包络线拟合问题、模态混叠等，这些问题导致HHT分析得到的物理信息不够理想。本文主要针对HHT过程中EMD分解和Hilbert变换所产生的端点效应问题进行研究。

1　端点效应分析

HHT在处理非平稳信号时存在以下两类端点效应问题：

第一类是EMD分解过程中出现的端点飞翼现象。在每一次的EMD筛分过程中，应根据信号的上、下极值点三次样条插值包络线来计算信号的局部平均值。由于信号不是无限长的，所以不能同时确定信号左右两端的局部极值点。如果直接在端点处弃值，将导致三次样条插值得到的上下包络线在信号两端附近扭曲失真(如图1)，产生端点飞翼。

图1　EMD分解端点飞翼示意图　　Fig.1　The end swing result of EMDdecomposition

针对EMD分解的端点飞翼现象，Huang等人最早提出利用原始信号进行延拓的思路。此后，许多文献提出用不同延拓的方法抑制端点飞翼：如极值点对称延拓[6]、多项式拟合延拓[7]、极值平移法[8]等均采用向原始数据添加极值点的思路；又如神经网络[9]、支持向量回归机[10-11]、最大相关波形[12]等方法则采用延长原信号数据的思路。

第二类是Hilbert变换过程中出现的端点发散现象[9]。对各个IMF分量进行Hilbert变换时，由于Hilbert变换是基于傅立叶变换实现的，而傅立叶变换会产生频谱泄露，出现严重的端点发散现象(如图2)，所得Hilbert谱无法真实地反映原始信号的特征。

图2　Hilbert变换端点发散示意图Fig.2　The end effects result of Hilbert transform

针对Hilbert变换过程中产生的端点发散现象，目前主要通过神经网络[9]、时变参数ARMA[13]和支持向量回归机[10]等数据延拓法先对IMF分量进行端点延拓，再进行Hilbert变换，并通过抛弃延拓的两端数据，使端点效应释放到原始数据的外端，保证原始数据范围内Hilbert变换的有效性。

上述HHT端点处理方法对特定信号的端点抑制效果较理想，但对于非线性、非平稳信号，由于其复杂与多变性，上述方法在计算精度和计算速度上难以同时达到最佳。本文将延长原信号数据与添加极值点相结合，提出在EMD分解过程中实施端部信号与极值点双重对称延拓的方法以抑制HHT的端点效应。与现有的极值点对称延拓等方法相比，这种新方法的端点效应抑制效果更佳，并能改善信号的整体正交性。与神经网络、时变参数ARMA和支持向量回归机的数据延拓等方法相比，本文算法更为简单，计算速度更快。并且，这种新方法不需要在Hilbert变换之前再对IMF分量进行端点延拓，可以在一次延拓中同时抑制EMD分解和Hilbert变换产生的端点发散。

2　双重对称延拓的端点效应处理

双重对称延拓首先进行端部数据对称延拓，按照原信号在端点附近的自然趋势延长原信号数据，把端点效应释放至原始信号外侧，从而使端点效应更缓慢影响到原信号内部。但是，如图3所示，如果仅进行端部数据延拓，三次样条插值所形成的包络曲线在端部依然存在端点发散，经过EMD分解，将影响到原始信号内包络曲线的准确性。因此，本文提出在端部数据延拓的基础上进一步进行极值点对称延拓，进一步控制包络曲线的走势，确保其准确性。

图3　端部数据对称延拓端点飞翼示意图Fig.3　The end flying result of end signal data symmetrical extension

2.1端部数据对称延拓

设原始信号为x(t)，对x(t)执行如下操作：

1)求x(t)的极大值indmax(1)，indmax(2)，…，indmax(n)和x(t)的极小值indmin(1)，indmin (2)，…，indmin (m)。其中，m=n，或丨m-n丨=1。

2)左端信号数据对称延拓：

当与左端点x(0)相邻的极值点为极大值时，比较x(0)与indmin(l)：

a.若x(0)>indmin(1)，以indmax(l)为对称点，将indmax(l)～indmax(3)之间的信号数据向左对称延拓(如图4(a)所示)。

b.若x(0)≤indmin(1)，以x(0)为对称点，将x(0)～indmin(2)之间的信号数据向左对称延拓(如图4(b)所示)。

当与左端点x(0)相邻的极值点为极小值时，比较x(0)与indmax(l)：

c.若x(0)

d.若x(0)≥indmax(1)，以x(0)为对称点，将x(0)～indmax(2)之间的信号数据向左对称延拓(如图4(d)所示)。

图4　左端数据对称延拓Fig.4　Left end data symmetrical extension

3)同理，对右端信号对称延拓，形成信号记为x1(t)。

2.2极值点对称延拓[6]

对经过端部数据对称延拓的信号x1(t)执行极值点对称延拓。比较信号的端点值与端部极值，当端点值大于相邻的第一个极大值或者是小于相邻的第一个极小值时，以端点作为对称点，当端点值在相邻极值中间时，以相邻极值点作为对称点，向信号外侧镜像添加两个极值点。

2.3按双重延拓后信号的包络线执行EMD筛分

1)确定信号x1(t)延拓的所有局部极值点，用三次样条曲线连接信号x1(t)所有的局部极大值和局部极小值，分别形成上包络线x1max(t)和下包络线x1min(t)。

2)求出x1max(t)和x1min(t)的平均值，记为m11(t)，将原数据序列x1(t)去掉该平均值后得到新数据序列h11(t)：

3) 判断h11(t)是否满足极值点个数与过零点个数相等或至多相差一个、局部极大值和局部极小值定义的包络均值为零这两个IMF的基本要求[1]，如果不满足则将h11(t)作为原始数据，重复上述1)～3)处理过程，直到新数据序列

满足IMF的条件，得到第一个IMF分量c1(t)：

4) 从x1(t)中分离出c1(t)，得到剩余序列r1(t)：

5) 将r1(t)作为一个新的原序列，按照以上1)～4)步骤，依次提取第2、3直至第n个固有模态函数cn(t)。当残量rn(t)成为一个单调函数或小于某一预定值时，分解结束。截取原始信号x(t)时域内IMF分量。

2.4EMD分解得到IMF分量直接进行Hilbert变换

信号x1(t)经EMD分解得到IMF分量ci(t)，不截取原始信号时域内的数据，保留延拓部分，直接进行Hilbert变换，将Hilbert变换的端点效应释放到原始数据外端：

构造解析信号

进一步取得瞬时频率

最后截取原始信号x(t)时域内瞬时频率。

3　仿真信号HHT分析

3.1EMD分解端点效应处理

已知简谐叠加信号x(t)=cos(2πt/50)+ 0.6cos(2πt/25) + 0.5cos(2πt/200)，t∈[1,300]。如图5所示。现按双重对称延拓方法进行HHT分析。

图5　简谐叠加信号Fig.5　Simple harmonic superimposed signal

图6是未经端点处理，直接以端点值作为极值点得到的IMF分量，其中虚线表示真实分量。该图显示，IMF1的端部信号明显收缩；随着频率减小，低频分量由于其时间尺度大，极值点之间的距离大，分量信号的内部更容易受到端部边缘效应的影响，导致EMD分解得到的低频分量失真。

图6　采用端点值作为极值点得到的IMF分量Fig.6　The IMF components using endpoint value as extreme point

图7为仅对该信号进行极值点对称延拓所获得的IMF分量。虽然端点效应在高频成份中得到了抑制，但在低频分量中，IMF2和IMF3右端点与真实值的误差达到30%和51%。

图7　采用极值点对称延拓得到的IMF分量Fig.7　The IMF components with extreme symmetrical extension

图8为仅对该信号进行端部数据对称延拓获得的IMF分量，该方法由于利用了更多的信号数据参与EMD分解，因此抑制效果略优于极值延拓的情况。但随着频率减小，低频分量的端点效应逐渐传播到信号内部，IMF3分量端点与真实值的最大误差仍达到30%。

图8　采用端部数据对称延拓得到的IMF分量Fig.8　The IMF components with signal end data symmetrical extension

图9为利用本文提出的双重对称延拓方法得到的IMF分量，与图8相比，该方法不仅在高频分量与真实分量吻合较好，低频分量也几乎接近真实信号。说明其更加有效地抑制了EMD分解的端点效应。

图9　采用双重对称延拓方法得到的IMF分量Fig.9　The IMF components with dual symmetrical extension

3.2EMD端点效应抑制效果评价指标

采用运算时间t、整体正交性指数IO[14]、信号分量相关系数ρ[15]等3个指标来进一步综合评价端点值作为极值点、极值点对称延拓、端点数据对称延拓、双重对称延拓等端点效应处理方法的抑制效果。

从表1可得，本文提出的双重对称延拓处理法的正交性指数IO值最小，端点值作为极值点的方法得到IO值最大，说明本文提出的双重延拓法比仅经端部数据延拓或极值点延拓的方法更有助于改善IMF分量的正交性。通过运算时间的比较可知虽然本文方法用时最长，但运行时间仅0.312 s，具有很高的计算效率。并且，本文方法的各个分量相关系数ρ最高，说明EMD分解结果与真实分量的误差最小。因此，经上述综合评价，双重对称延拓法是一种端点效应抑制良好、算法简单、计算高效、IMF分量之间正交性程度较高的改进方法。

表1不同延拓方法得到的端点效应抑制效果评价指标

Tab.1The evaluation indexes of end effects reduction by different extension methods

评价指标端点值作为极值点极值点对称延拓端部数据对称延拓双重对称延拓IO值0.05760.03780.03360.0199用时/s0.28000.29700.29200.3120ρ10.93950.99510.99670.9982ρ20.95760.99160.99810.9993ρ30.88370.98100.93890.9974

3.3Hilbert变换的端点效应处理

将上述仿真信号x(t)采用端点值作为极值点方法得到瞬时频率如图10所示，端点发散效应清晰可见。

图10　采用端点值作为极值点得到的瞬时频率Fig.10　The instantaneous frequency with endpoint value as extreme point

利用极值点对称延拓后得到的瞬时频率见图11，边缘效应虽有所改善，但仍然没有完全得以抑制。

图11　采用极值点对称延拓得到的瞬时频率Fig.11　The instantaneous frequency with extreme symmetrical extension

图12是采用端部数据对称延拓方法得到的瞬时频率，其端点问题释放到原始数据的外端，其端点发散有显著的抑制，仅在瞬时频率为0.005 Hz的曲线右端出现发散，发散部分与真实值的最大误差达到45%。可见，数据延拓的思路对于抑制EMD分解和Hilbert变换两阶段的端点效应均有效，但在低频分量中，仍可能因延拓部分包络失真影响到原始信号端部的准确性。

图12　采用端部数据对称延拓得到的瞬时频率Fig.12　The instantaneous frequency with signal ends data symmetrical extension

图13是采用本文提出的双重对称延拓方法得到的瞬时频率图，从图中可清楚看出3个分量频率分别对应仿真信号的3个频率成份0.04、0.02、0.005 Hz。相比，瞬时频率为0.005 Hz的曲线右端更加平稳，端点效应抑制得更加彻底。

图13　采用双重对称延拓得到的瞬时频率Fig.13　The instantaneous frequency with dual symmetrical extension

4　工程应用

现通过结构地震响应的分析进一步验证本文方法的有效性。已知单自由结构体系质量为107kg，刚度为5×107N/m，结构自振频率为0.355 9 Hz，结构阻尼比为0.05。输入荷载为cos 2πt，采样时间间隔为0.01 s，采样的个数为1 000。经弹性时程分析，该结构绝对加速度响应如图14所示。

图14　结构绝对加速度响应Fig.14　Structural absolute acceleration response

信号经过双重对称延拓处理后的EMD分解结果如图15所示。由图可知绝对加速度响应由高频至低频依次分解得到3阶IMF分量和一个剩余量，说明EMD分解中分辨率是自适应的。各个分量两端未出现大幅度的摆动现象，这说明提出的双重延拓方法能较好的抑制EMD端点飞翼。

图15　双重对称延拓处理的EMD结果Fig.15　EMD result with dual symmetrical extension

将保留延拓数据的IMF分量进行Hilbert变换得到的瞬时频率如图16所示，从图中可以清晰看出前三阶分量的频率分别约为1、0.36、0.2 Hz，同时频率信号端部也无明显的端点发散现象，说明双重对称延拓对于实际应用来说也是有效的。EMD分解出来的IMF分量可以得到原始信号无法得到一些物理信息。图15中的前两个分量描述了体系动力响应的本征振动模态：第一个分量描述的是体系动力响应中与输入荷载有关的稳态振动部分；第二个分量描述的是体系动力响应中与结构自振特性有关的伴生自由振动部分。

图16　前三阶IMF分量对应的瞬时频率Fig.16　The instantaneous frequency of the first three order IMFs

表2给出了4种分解方法对于结构绝对加速度响应处理后的两个性能评价结果，其中双重对称延拓的IO值最小，EMD分解的各个分量正交性最好。同时运行的时间只有0.3 s，运行速度快。

表2绝对加速度响应信号不同延拓对比

Tab.2The comparison of extension among absolute acceleration response(signals)

评价指标端点值作为极值点极值点对称延拓端部数据对称延拓双重对称延拓IO值0.06810.06330.06520.0520用时/s0.22000.28000.26000.3000

5　结论

以双重对称延拓的方法较端点值作为极值点、极值点对称延拓、端部数据对称延拓等方法能更好地抑制EMD分解带来的端点飞翼问题，同时改善了信号的整体正交性，使Hilbert变换的端点效应释放到原始信号的外端，从而抑制了Hilbert变换的端点发散效应。与神经网络、时变参数ARMA和支持向量回归机的数据延拓等方法相比，本文提出的方法更为简单，计算速度更快，可以在一次延拓中同时抑制EMD分解和Hilbert变换产生的端点发散。同时该方法在实际信号中也能有效地抑制端点效应问题。

参考文献：

[1] Huang N E, Shen Z, Long S R.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of Royal Society of London,1998,454: 903-995.

[2] 郑近德，程军圣.改进的希尔伯特-黄变换及其在滚动轴承故障诊断中的应用[J].机械工程学报,2015(1):138-145.

[3] 熊炘.基于希尔伯特-黄变换的故障转子振动模式分析方法研究[D].杭州：浙江大学,2012.

[4] 吴祚菊，张建经，王志佳，等.改进HHT算法在地震动特性分析中的应用研究[J].地震工程与工程振动，2015(2):86-93.

[5] 吴琛，周瑞忠.Hilbert-Huang变换在提取地震信号动力特性中的应用[J].地震工程与工程振动,2006(5):41-46.

[6] Rilling G, Flandrin P, Goncalves P.On empirical mode decomposition and its algorithm[C]//Proceeding of IEEE-EURASIP Workshop on Nonlinear Signal and Image (NSIPL2003).Trieste: Institute of Electrical and Electronics Engineers,2003:8-11.

[7] 刘慧婷，张旻，程家兴.基于多项式拟合算法的EMD端点问题的处理[J].计算机工程与应用,2004,40(16):84-86,100.

[8] 王学敏，黄方林.EMD端点效应抑制的一种实用方法[J].振动、测试与诊断，2012(3):493-497.

[9] 邓拥军，王伟，钱成春，等.EMD方法及Hilbert变换中边界问题的处理[J].科学通报，2001(3):257-263.

[10] 程军圣，于德介，杨宇.基于支持矢量回归机的Hilbert-Huang变换端点效应问题的处理方法[J].机械工程学报，2006(4):23-31.

[11] 旷欢，王如龙，张锦，等.基于SVM的EMD端点效应抑制方法研究[J].计算机工程与应用，2015(11):196-200.

[12] 高强，段晨东，赵艳青，等.基于最大相关波形延拓的经验模式分解端点效应抑制方法[J].振动与冲击，2013(2):62-66.

[13] 程军圣，于德介，杨宇.Hilbert-Huang变换端点效应问题的探讨[J].振动与冲击，2005(6):40-42.

[14] 肖韵.Hilbert-Huang变换端点问题处理方法的分析与研究[D].成都：西南交通大学,2012.

[15] 王秋生，段丹辉.经验模态分解的边界效应处理技术[J].计算机测量与控制，2006(12):1673-1675.

(责任编辑： 陈雯)

A dual symmetrical extension method for end effects reduction of Hilbert-Huang transform

Xiang Hong1,2, Wu Chen1,2, Du Xipeng1,2

(1.College of Civil Engineering, Fujian University of Technology, Fuzhou 350118, China;2.Fujian Provincial Key Laboratory of Advanced Technology and Informationization in Civil Engineering, Fuzhou 350118, China)

Abstract：End effects reduce the precision of empirical mode decomposition (EMD) greatly．A dual symmetrical extension method was proposed to decrease the end effects of EMD.Signal end data symmetrical extension was regarded as the first extension, whereas extreme symmetrical extension was regarded as the second extension.The proposed method can also decrease the end effects of Hilbert transform.The harmonic superimposed signal and absolute acceleration response were applied to test the performance of the proposed method.The results show that the proposed method can reduce the end effects of HHT and can improve its precision significantly.

Key words：Hilbert-Huang transform; empirical mode decomposition; end effect; dual symmetrical extension method; evaluation index

doi:10.3969/j.issn.1672-4348.2016.04.004

收稿日期:2016-05-31

基金项目：国家自然科学基金项目(51108083); 福建省自然科学基金项目(2015J01631)

第一作者简介：项洪(1991-)，男，福建龙岩人，硕士研究生，主要从事结构抗震理论与工程应用研究.

中图分类号：TN911.72

文献标志码：A

文章编号：1672-4348(2016)04-0324-08