EMD 包络线拟合算法改进及在泡沫尺寸趋势提取中的应用

卢 青， 王雅琳 ， 彭 凯， 薛永飞， 谢永芳

(中南大学 信息科学与工程学院，湖南 长沙410083)

经验模态分解法(Empirical Mode Decomposition，EMD)在非平稳及非线性信号的处理及分析上具有明显优势，被广泛应用于海洋内波探测［1］及温度预测［2］、大气时间序列预测［3］、天体观测［4］、地震记录［5］、机械故障诊断［6］等方面。包络线拟合是EMD 算法中非常关键的一步。目前，常用的包络线拟合方法是三次样条插值法［7-10］。但现有三次样条插值拟合包络线的方法存在着插值点选取不当、端点发散［11］，拟合松散［7］等问题，致使拟合误差的出现，该拟合误差会随着筛分过程逐渐向内，使得分解的整个数据被严重污染，导致分解结果失真。

目前，在包络线拟合方面的研究主要是针对端点发散问题，在端点效应抑制方面国内外学者提出了很多种方法，具体可归为3 类［12］:波形延拓法、数据预测延拓法和极值延拓法。波形匹配［13］、波形镜像等方法属于波形延拓法，神经网络预测［14］、ARMA 模型预测［15］、时间序列预测等方法是数据预测延拓，多项式拟合、极值对称延拓法［11］等方法属于极值延拓法。其中，常用的端点效应抑制方法为波形匹配、神经网络预测、ARMA 模型预测、极值对称延拓。然而波形匹配法计算量大，占用内存;神经网络预测方法的学习过程时间较长，大大影响了EMD 的分解速度，不适合实时信号处理;ARMA 模型用于处理非线性非平稳的信号时参数估计比较复杂，且分析效果受限;极值对称延拓是镜像延拓的延续与提升，方法简单，但处理非周期信号时误差很大;其余各种方法针对非线性非平稳信号延拓结果因人而异，很难保证广泛的适用性。因此需要提出一种简单、容易实现、适用性强、有效的方法。

文中通过加权最小二乘多项式拟合方法进行极值延拓，根据延拓的实际情况对得到的端点值进行合理性判断并加以取舍，然后作为插值拟合曲线的端点;同时结合包络线拟合中存在的另外两个问题，提出了EMD 中包络线拟合的改进方法。为评价改进方法的包络线拟合效果，文中定义了包络线与信号的方差和以及曲线长度协同作为评判标准。

铝土矿泡沫浮选中不同槽间的泡沫尺寸变化趋势存在关联规则，要挖掘相关规则，需先提取每个槽的泡沫尺寸的变化趋势，从而获取泡沫尺寸变化的趋势图，为关联性分析打下基础。文中将所提包络线拟合改进方法应用于粗选槽泡沫尺寸的趋势提取，获得了较好的粗选槽泡沫尺寸变化趋势，验证了所提方法的有效性。

1 EMD 算法及其存在的问题

1.1 EMD 算法简介

美籍华人HUANG 等于1998 年提出了一种适合于分析非线性、非平稳信号序列的经验模态分解法EMD。该方法可以将复杂信号分解为有限数量的含原信号不同时间尺度的局部特征信号的本征模函数(Intrinsic Mode Function，IMF)［3］，进而根据不同频率分量的特性进行滤波或预测。

EMD 算法步骤如下:

1)寻找信号X(t)的所有极值点。

2)分别以所有极大值点、极小值点为插值点进行三次样条插值，得到信号X(t)的上、下包络线。

3)求出上下包络线的平均值m(t)，用以表示信号的平均值曲线。

4)X(t)－m(t)= h1，完成一次筛分过程。第2次筛分过程用h1代替X(t)作为分解信号，重复k次，直到h1k满足IMF 条件:①函数在整个时间范围内，局部极值点和过零点的数目必须相等，或最多相差一个;②在任意时刻，局部最大值的包络(上包络线)和局部最小值的包络(下包络线)平均为零。

5)记c1= h1k，做运算X(t)－ c1= r1，再将r1作为分解数据，重复筛分过程，得到c2，r2，…，cn，rn，直到rn或cn小于给定值或rn成为单调函数不可再筛分出IMF。则分解完成，有

其中，ci为第i 个IMF 分量，代表了原始信号X(t)中不同特征时间尺度的信号分量;rn为残余函数，代表信号的平均趋势。

1.2 存在问题及原因

EMD 的关键步骤之一是用三次样条函数拟合得到信号的上、下包络，它将直接影响分解得到的所有IMF 分量的精度。然而在进行信号的包络线拟合时，由于信号在端点位置往往不是极值点，使得拟合的包络线在边界附近出现严重失真，这种现象被称为边界效应或端点效应。实际应用中信号中存在的噪声及偏差不可能完全去除，可能导致极值点数量增加;同时由于边界效应产生的误差及极值点数量的增加将引起拟合误差，且随着筛分过程的不断反复逐渐传播到信号的内部并污染整个信号，最终导致得到不可靠的IMF 分量和残余分量。

例如，分段信号

图1 是用三次样条插值方法在时间段［0，20］内拟合信噪比为70 的信号y0(t)的波形图。

图1 三次样条插值法拟合y0(t)的包络线Fig.1 Envelope extraction of y0(t)using CSI

包络线应当完整地包络信号，且是光滑的曲线，上包络线的所有值不小于对应点的信号值，下包络线的所有值不大于对应点的信号值。由图1 可以看出，三次样条插值法对于非平稳信号的拟合效果不佳。究其原因:

1)所采用的信号不是平滑信号。实际应用中由于采集设备误差、记录误差及现场环境影响等信号数据本身就不完全平滑，且信号采集过程中无法避免噪声信号的干扰，滤波和去噪在一定程度上可以净化信号，但提取包络线时信号中必然还残留一些噪声，容易出现极值点冗余现象，若把这些极值点全部作为插值点，会出现包络线贴合于原信号曲线的现象，导致曲线拟合失真。

2)实际应用中采集到的信号可能会出现某一小段时间内幅值不变的情况，由极值点的定义可知这一时间段内没有极值点，若依旧只是把极值点作为插值点就会出现一段包络线的空缺。

3)端点无约束，包络线在端点处出现发散现象。

4)相邻的两个插值点之间的插值曲线斜率出现过零点，即拟合曲线不单调，部分包络值在临近极值之外，出现了“凸峰”或“凹谷”现象，导致包络线拟合的不“紧”。

2 算法改进及效果评判

针对1.2 中提到的包络线拟合中插值点的选取、端点效应的抑制以及插值曲线拟合松散的问题，经综合考虑，提出有效的解决办法对EMD 中包络线的拟合予以改进。

以同时满足包络线与信号的方差和尽可能小且包络线长度尽可能小，作为评判包络线拟合效果的标准。包络线与信号的方差和小说明包络线与信号曲线贴合度好，无过冲现象及端点发散现象，但包络线与信号的方差和太小又说明包络线和信号有重叠趋势，包络线拟合效果不好，因此以包络线与信号的方差和单独作为评价标准有一定合理性却并不科学。包络线长度值小说明无过冲、欠冲、端点发散现象，但包络线长度最短，即以直线连接各插值点时，不满足包络线光滑的要求，因此包络线长度最短单独作为评价标准也有一定合理性却并不科学。所以，以两项指标加权和e 作为评判包络线拟合效果的标准。权值选择的原则是使两个标准计算值的数量级相同，即es+ex与| ys| +| yx| 的数量级相同，选取办法是将数量级高的一个标准的数值降级。

记上包络线为ys(t)，下包络线为yx(t)，上、下包络线与信号的方差和分别为es和ex，上下包络线的长度分别为| ys| 和| yx|，则由方差和及曲线长度的定义可得

式中:i 为第i 个上包络线插值点;p 为上包络线插值点个数;ysi为第i 与i + 1 个插值点间的插值函数;y'si为第i 与i +1 个插值点间插值函数的一阶导数，i = 0 时取上包络线左端点，i = p 时i +1 = p +1 取上包络线右端点;j 为第j 个下包络线插值点;q 为下包络线插值点个数;yxj为第j 与j +1 个插值点间的插值函数;y'xj为第j 与j +1 个插值点间插值函数的一阶导数，j = 0 时取下包络线左端点，j = q 时j +1 =q +1 取下包络线的右端点。

对于y0(t)，鉴于包络线与信号的方差和是包络线长度数量级的100 倍，所以取

图1 中三次样条插值法直接拟合曲线y0(t)的误差见表1。

表1 三次样条插值法拟合包络线的误差Tab.1 Error of the envelope extraction using CSI

2.1 基于特殊点的三次样条插值

首先，针对小段值不变的平稳数据导致的无极值点问题，提出一种基于特殊点的三次样条插值方法。特殊点既包括传统意义上的极值点又包括平稳数据段的所有点。具体选取方法为:上包络线插值点x(ti)满足

下包络线插值点x(tj)满足

其中，M，N 分别为上、下包络线插值点的个数。

其次，针对信号误差引起的小段时间内极值点数量增多且极值波动很小导致拟合失真的问题，选择插值点时，在寻找到特殊点后增加一步筛选操作。以下包络线的拟合为例阐述筛选原理:将所有极小值点组成新的时间序列，按照一定的时间间隔Δ 对x1(t)进行分段，对每一段的所有时间数据求均值，该时间均值处的点作为插值点进行包络线拟合。具体步骤如下:

1)由时间序列的时间特性可知，后一时刻的时间值大于前一时刻的时间值，即ti+1＞ti＞ti－1。将时间序列x1(t)在ti－1与ti之间分段，分界点产生规则是:①ti与前一点ti－1的时间间隔大于一定的阈值σ，即Δ = ti－ti－1＞σ，σ 由时间序列本身的波动周期(周期序列)及序列长度决定。②若x1(ti)＞x1(ti+1)，则必须成立;否则，若x1(ti)＜x1(ti+1)，则必须成立。阈值R由序列本身决定。

2)对每一段数据的时间坐标求算数平均¯t，用该均值处的点(¯t，x1(¯t))近似作为该段的唯一插值点。改进后基于特殊点的包络线拟合如图2 所示，拟合误差见表2。

图2 基于特殊点的包络线拟合Fig.2 Envelope extraction based on the special points

表2 基于特殊点的包络线拟合误差Tab.2 Error of the envelope extraction based on the special points

通过图2 及表2 与图1 及表1 的对比可发现，3个信号拟合包络线的长度及包络线与信号的方差和都有明显减小，误差e 均大幅减小，从图2 中也可以看出拟合效果有明显改善。

2.2 端点效应的抑制

端点发散是因为端点无约束，所以解决办法是增加约束。若获取的信号量足够多，采样时在有用信号的前后分别多采集一段数据，以增加极值点的个数作为约束，这样得到的包络线在端点处不会严重扭曲。若已知曲线方程，可以在信号前后再计算出几个极值点作为插值点进行包络线拟合。但是在没有多余数据又不知道曲线方程的前提下，在拟合包络线时，需要对已知插值数据进行延拓，以在数据两端获得附加的极大值和极小值作为约束抑制端点效应。

加权最小二乘多项式拟合是对所有的点进行整体拟合，而三次样条插值是在相邻两个插值点间建立三次多项式进行拟合的分段拟合，端点处为自由端点，所以加权最小二乘多项式拟合得到的端点比自由端点精确得多。为此，提出加权最小二乘多项式拟合实现EMD 包络线拟合中端点发散的抑制。

以下包络线左端点的确定为例说明采用加权最小二乘多项式拟合确定约束条件的方法。加权最小二乘原理是找到一个函数

可以最佳地拟合于诸观测点

即应当满足

此处选取拟合多项式的最高次数为3 次，尽管拟合次数越高，拟合精度越高，拟合效果越好，但4次拟合计算量急剧增大，且易出现放大微小误差的问题。且3 次拟合与2 次拟合相比，计算量相差不多，拟合精度却会有所提升。相应地，插值点至少选取4 个。但当特殊点数量不足4 个时，选取全部的特殊点，并相应地可降低拟合多项式的次数。此时

需满足

权值选取原则是误差项方差越大则权值越小，权值的最佳取值是方差和的倒数。但是由于方差和不可知，习惯选取方法是取横坐标x 的倒数。此处，根据实际情况，x 应理解为观测点横坐标xi与左端点横坐标x0的距离，即

若增加的下包络线端点值大于临近的极大值点，显然不符合实际，舍弃该端点不用，以靠近端点的第一个下包络特殊点和第一个上包络特殊点的均值作为端点值。虽然这样取得的端点值存在误差，但由于本意并不是获得精确的端点值，只要可以有效抑制端点发散问题就是可取的。同理，对上包络线添加端点，如图3 所示的绘出包络线，拟合误差见表3。

图3 抑制端点效应后的包络线Fig.3 Extraction with the restraining end effect

表3 抑制端点效应后的拟合误差Tab.3 Error of extraction with the restraining end effect

由图3 及表3 与图2 及表2 相比较可以看出，发现拟合误差e 减小，曲线拟合效果有一定的改善。

2.3 保形分段三次Hermite 插值

由图3 可以看出，在抑制端点效应之后三次样条插值法拟合出的包络线依然存在与信号曲线有交点(下包络线部分点取值大于信号值或上包络线部分取值小于信号值)，且相邻插值点间的下包络线有过冲现象，拟合的包络线松散。原因是相邻两插值点间存在单调性不一致，包络线缺少约束。

针对上述问题，首先要增加一组对包络线拟合的约束条件:上包络线值始终不小于相同坐标处的信号值;同时下包络线的所有值都不大于相同坐标处的信号值。即对于所有时间点t，有ys(t)≥y(t)≥yx(t)。然后采用分段三次埃尔米特多项式插值［16］(Cubic HermiteInterpolation，CHI)进行包络线拟合，可以有效改善这一缺陷。

在区间［a，b］上构造三次埃尔米特插值时，首先取a = x0＜x1＜x2＜… ＜xn= b，记xi处函数值yi，导数值y'i;然后在每个子区间［xi，xi+1］，(i =0，1，…，n － 1)上构造三次插值多项式Hi(x)，Hi(x)需满足Hi(xi)= yi，H'i(xi)= y'i。但一阶导数y'i通常不可知，所以y'i的确定就成为CHI 的关键。

具体步骤为:

1)借鉴Lagrange 插值基函数法，引入4 个三次多项式α0(x)，α1(x)，β0(x)，β1(x)作为基函数，令Hi(x) = yiα0(x) + yi+1α1(x) + y'iβ0(x) +y'i+1β1(x)，并求解基函数。

2)于是，x ∈［xi，xi+1］时，两点Hermite 插值

所以，［a，b］上的插值函数为

3)一阶导数的确定:(xi，yi)(i = 1，2，…，n －1)前后子区间(xi－1，xi)，(xi，xi+1)的斜率分别为

所以(xi，yi)点处斜率di为

其中

CHI 插值具有在每个子区间单调的特性，且公式简单、计算量小、稳定性好［17］。基于CHI 插值的包络线拟合如图4 所示，拟合误差见表4。

图4 CHI 拟合包络线Fig.4 Envelope extraction by using CHI

表4 CHI 拟合后的误差Tab.4 Error of extraction using CHI

与图3 及表3 相比，误差e 减小，同时可以看出包络线在两插值点间拟合曲线的过冲、欠冲，即拟合松散问题得到改善。由此证明了文中所提出的改进方法有效。

2.4 仿真验证

图5(a)、(b)分别为信噪比为30 的非线性非平稳信号

直接用三次样条插值法拟合及改进方法拟合的包络线，拟合误差见表5。

图6(a)、(b)分别为信噪比为30 的信号y2(t)= sin(8πt)+5sin(0.4πt)+2sin(2πt +π/3)，直接用三次样条插值法拟合及改进方法拟合的包络线，拟合误差见表6。

图5 y1(t)拟合包络线Fig.5 Envelope extraction of y1(t)

表5 y1(t)拟合误差对比Tab.5 Comparison of extraction error of y1(t)

图6 y2(t)拟合包络线Fig.6 Envelope extraction of y2(t)

表5 和表6 均取

通过上述两例的对比可以发现，e 减小，同时可以看出改进后的包络线拟合方法拟合出的包络线包的更紧密，端点发散现象得到抑制。由此验证了所提改进方法在包络线拟合上的有效性。

表6 y2(t)拟合误差对比Tab.6 Comparison of extraction error of y2(t)

3 改进的包络线提取方法在泡沫尺寸变化趋势提取上的应用

典型的铝土矿泡沫浮选操作流程，包括粗选、精选、粗扫和精扫作业。实际生产中，操作工人通过观察泡沫进行生产操作。工业上通过分布机器视觉系统在每个工序获取泡沫图像并从中提取出颜色、尺寸、纹理等信息。鉴于各工序泡沫存在关联信息，对不同工序之间的泡沫信息进行关联性分析得到关联规则，可以用于指导生产，对泡沫尺寸的包络线分析可以获取泡沫尺寸变化的规律，便于进行关联性分析。

选取基于分布机器视觉的铝土矿泡沫浮选过程在粗选槽获取的图像中提取的一班泡沫尺寸数据，时间间隔为5 min，剔除异常值后作为原始数据。针对这组数据，直接用三次样条插值法拟合包络线与用改进后的方法拟合的包络线如图7 所示，拟合误差见表7。

对比图7(a)和图7(b)及表7 中所列拟合误差，可以发现包络线与信号的方差和以及曲线长度同步减小，误差e 减小;同时可以看出改进后的包络线拟合方法拟合出的包络线包的更紧密，不存在插值点冗余及端点发散问题。再次验证了所提改进方法在包络线拟合上的有效性。

表7 包络线拟合方法改进前后的拟合误差Tab.7 Extraction error before and after improvement

4 结 语

文中分析了EMD 算法中基于三次样条插值的包络线拟合方法，列出了存在的问题并分析了原因，提出了基于特殊点的EMD 包络线拟合改进算法。在极值点的基础上增加了转折点与值不变的平稳点组成特殊点，经过分段筛选删除冗余插值点。针对端点发散问题，通过加权最小二乘多项式拟合对所有极值点进行多项式拟合计算出端点值。三次埃尔米特插值基于其在每个子区间单调的特性及对上下包络线数值的约束可以有效改善三次样条插值包络线拟合中存在的包络线松散且与信号有交叉的情况。

在有针对性的提出相应改进方法的同时，通过定义的拟合优良性评判标准，及逐步绘出改进后的包络图，构成清晰明确的客观数据及主观视觉的对比，表明了所提改进方法在包络线拟合上的良好效果。最后将改进方法应用于铝土矿泡沫浮选粗选槽获取的泡沫尺寸数据，提取其变化趋势，再次证明了所提方法对包络线拟合的有效性。

［1］宋海斌，拜阳，董崇志，等. 南海东北部内波特征——经验模态分解方法应用初探［J］. 地球物理学报，2010，53(2):393-400.

SONG Haibin，BAI Yang，DONG Chongzhi，et al.A preliminary study of application of empirical mode decomposition method in understanding the features of internal waves in the northeastern south china sea［J］.Chinese Journal of Geophysics，2010，53(2):393-400.(in Chinese)

［2］Bhuiyan S M A，Carey S，Khan J F.Trend analysis of sea surface temperature and sea level pressure employing a pseudo-EMD method［C］//2013 Proceedings of IEEE Southeastcon.Jacksonville，FL:IEEE，2013:1-4.

［3］玄兆燕，杨公训.经验模态分解法在大气时间序列预测中的应用［J］.自动化学报，2008，34(1):97-101.

XUAN Zhaoyan，YANG Gongxun.Application of EMD in the atmosphere time series prediction［J］.Acta Automatic Sinica，2008，34(1):97-101.(in Chinese)

［4］XU T，WU J，WU Z S，et al.Long-term sunspot number prediction based on EMD analysis and AR model［J］.Chinese Journal of Astronomy and Astrophysics，2008，8(3):337-342.

［5］钱昌松，刘代志，刘志刚，等. 基于递归高通滤波的经验模态分解及其在地震信号分析中的应用［J］. 地球物理学报，2010，53(5):1215-1225.

QIAN Changsong，LIU Daizhi，LIU Zhigang，et al.EMD based on recursive high-pass filter and its application on seismic signal analysis［J］.Chinese Journal of Geophysics，2010，53(5):1215-1216，1224-1225.(in Chinese)

［6］程军圣，于德介，杨宇.基于内禀模态奇异值分解和支持向量机的故障诊断方法［J］.自动化学报，2006，32(3):475-480.

CHENG Junsheng，YU Dejie，YANG Yu. Fault diagnosis approach based on intrinsic mode singular valuedecomposition and support vector machines［J］.Acta Automatic Sinica，2006，32(3):475-480.(in Chinese)

［7］朱伟芳，赵鹤鸣，陈小平.一种最小长度约束的EMD 包络拟合方法［J］.电子学报，2012，40(9):1909-1912.

ZHU Weifang，ZHAO Heming，CHEN Xiaoping. Aleast length constrained envelope approach for EMD［J］. Acta Electronica Sinica，2012，40(9):1909-1912.(in Chinese)

［8］张勇，陈滨.黄变换的穿越筛分法［J］.电子科技大学学报，2007，36(1):8-10.

ZHANG Yong，CHEN Bin.Crossing sifting method based on huangtransform［J］.Journal of University of Electronic Science and Technology of China，2007，36(1):8-10.(in Chinese)

［9］JIA X F，AN H Q，ZHANG S G.Cubic spline interpolation method for theenvelope tracking of middle and low frequency voltage flicker［J］.Thermal Power and Electrical Engineering III，2014(960/961):704-709.

［10］QIU Mianhao，LIU Jing，CONG Hua.Research of direct sifting EMD based on cubic B-spline interpolation curve and its application in processing mechanical vibrating signals［J］.Journal of Academy of Armored Force Engineering，2007，21(3):29-33.

［11］GUO Yanwei，ZHOU Tao. The new method based on improvement EMD end effect for power quality analysis［J］. Electrical Engineering，2013，12:13-16.

［12］曹端超，康建设，赵劲松，等.EMD 端点效应抑制方法仿真比较与实例分析研究［J］.机械传动，2013，37(3):83-87.

CAO Duanchao，KANG Jianshe，ZHAO Jingsong，et al. Research on the comparison with methods of restraining endingeffect of EMD and its application in fault diagnosis［J］. Mechanical Transmission，2013，37(3):83-87.(in Chinese)

［13］李合龙，王龙，李明建，等.基于平均包络线匹配算法的EMD 端点效应分析及在股价趋势分解中的应用［J］.系统工程理论与实践，2013，33(8):2072-2079.

LI Helong，WANG Long，LI Mingjian，et al.The end effect of EMD based on matching mean envelope and itsapplications in trend decomposition of stock price［J］.Systems Engineering-Theory and Practice，2013，33(8):2072-2079.(in Chinese)

［14］孟宗，顾海燕，李姗姗.基于神经网络集成的B 样条经验模态分解端点效应抑制方法［J］.机械工程学报，2013，49(9):106-112.

MENG Zong，GU Haiyan，LI Shanshan. Restraining method for end effect of B-spline empirical modedecomposition based on neural network ensemble［J］.Journal of Mechanical Engineering，2013，49(9):106-112.(in Chinese)

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

CHENG Junsheng，YU Dejie，YANG Yu.Discussion of the end effects in Hilbert-Huang transform［J］.Journal of Vibration and Shock，2005，24(6):40-42，47.(in Chinese)

［16］邓林峰，赵荣珍.基于CHI-LMD 方法的转子振动信号分析［J］.振动与冲击，2014，35(15):58-64.

DENG Linfeng，ZHAO Rongzhen. Rotor vibration signal analysis based on the CHI-LMD method［J］. Journal of Vibration and Shock，2014，35(15):58-64.(in Chinese)

［17］ZHU W，ZHAO H，CHEN X. Improving empirical mode decomposition with an optimized piecewise cubic hermiteinterpolation method［C］//2012 International Conference on Systems and Informatics (ICSAI). Yantai:IEEE，2012:1698-1701.