二维经验模式分解算法的改进研究

刘 兵，李 聪，向磊磊，孙玉秋 （长江大学信息与数学学院，湖北荆州434023）

Norden Huang于1988年提出一种新的用于分析非线性和非平稳的信号处理方法——经验模式分解 （Empirical mode decomposition，EMD）［1］，这种方法的关键在于：任何复杂的数据都可以分解为一系列有限且是少量的固有模式函数 （Intrinc mode function，IMF），其中IMF满足以下2个条件：①在整个数据序列中，极值点的数量与过零点的数量必须相差不超过1个；②在任一时间点上，信号局部极大值确定的上包络线和局部极小值确定的下包络线的均值为零，即满足传统平稳高斯过程。这种分解是自适应的，因此具有很高的效率，将这种分解方式应用于信号分解时能获得比较好的IMF，并且对非线性和非平稳的信号进行处理时，能够取得小波分析等分析方法难以取得的效果。

EMD方法在处理一维信号时，对极值的选取一般只需与左右像素比较，在拟合时可以得到足够的拟合点。但是在二维图像中进行选取极值的过程，一般需要比较每个像素的3邻域，5邻域或8邻域，然后在拟合时对每列或每行进行拟合，可能会出现某行或某列的极值点低于2个，这样算法就不具有普遍性。应用形态学的算法可以很好的解决此类问题，并且所需的存储空间远低于原算法。在信号中，对极值点的拟合一般采用三次样条拟合，这种拟合方式会出现拟合过冲的情况，且在端点处也有可能会发生端点飞翼。可是这些变异点在信号处理中不是很多，所以没有引起研究者的重视。但是在图像处理时，这些变异点会随着计算次数的增加来影响整个图像［2］。针对这些问题，笔者提出1种EMD的改进算法：先运用形态学算子［3］提取出图像的极大值和极小值点，对极值点采用Delaunay三角拟合［4-5］算法进行曲面拟合，得到上下包络面。

1 算法改进

1．1 EMD算法

一维EMD算法原理如图1所示。把二维EMD算法［6］运用到图像处理方面，将一幅图像分解为若干固有模式函数和一个余量函数的集合，和一维情况类似，二维EMD定义如下：

1）把原始数据X作为待处理数据，确定该图像的所有局部极值点 （包括极大值点和极小值点）；

2）对所有极大值和极小值分别进行曲面拟合得到包络曲面emax，emin；

4）计算余量h＝X-eave；

5）判断h是否是一个固有模式分量 （IMF），若不是，将h返回到第一步。

经过第5）步的判定，假如条件满足，就可以得到一个固有模式函数，将原始图像减去所求得的IMF1，即得到一个余量数据r，将r作为原始数据再重复前面5个步骤，得到IMF2，IMF3，…，IMFn以及一个余量R。余量R满足预先设定的停止准则后即可停止，最后剩下原始数据的余项R。

将这些固有模式函数和余量相加，就可以实现对信号的重构。

图1 一维EMD算法流程

1．2 利用形态学提取极值点

用原始的方式进行极值点的提取，对于边界上的灰度值要单独处理，耗费大量时间且寻找的极值点在某一行 （列）只有一个或零个，对后续操作有难以克服的影响。因此，考虑用形态学的方法来求极值。

形态学有以下4种基本运算：膨胀，腐蚀，开运算，闭运算，用形态学求极值主要运用的是膨胀和腐蚀运算。所谓在灰度级图像上的膨胀操作就是指输入图像在结构元素 （本身可以看作是一个子图像函数）下的一种灰度扩张。原图像滑过结构元素时，图像有部分与结构元素重合，此时对应点的灰度值根据结构元素的取值进行取最大值处理。处理过后的原图像主要取决于结构元素的值和形状，如果所有结构元素的值为正，则输出图像会趋于比输入图像更亮并且暗的部分全部减少或者被消除掉。同样腐蚀操作是指输入图像在结构元素下的一种灰度降低。与膨胀运算相比，仅在对应点灰度值取值时，所采取的运算方式不一样，腐蚀操作是根据结构元素采取最小值操作。

选取一个3阶元素全为1的矩阵作为结构元素，运用形态学的膨胀 （腐蚀）算子，分别得到的是每一个点在其3邻域内的最大值 （最小值）为这一点的值。然后与原图像进行比较，具有相同数值的点就是原图像的极大 （小）值点。

1．3 利用Delaunay三角插值拟合

一维的插值算法在二维图像中直接应用时会产生较大的误差，因此，采用Delaunay三角插值算法，直接对极值点进行插值拟合，这样可以极大地减小图像因为拟合而引起的误差。

二维点集三角剖分是指将二维平面上的点集用不相交的直线段连接起来，使得所形成的凸包内每一个区域都是三角形。Delaunay三角插值算法是先用Voronoi多边形将离散的极值点分开，使每一个极值点属于一个Voronoi多边形，而连接3个共点的Voronoi多边形内的极值点则形成一个Delaunay三角形，所有这些Delaunay三角形的集合构成Delaunay三角剖分。在每个三角形剖分内进行三次样条插值，就得到了点集的拟合面。

Delaunay准则实现的条件是凸包内每一个三角形外接圆中不包含点集中的其他任何点，这使得每个三角形都尽可能接近于等边三角形，避免产生狭长的三角形也就是使各离散点对整个三角形有限元网格的影响仅限于局部。

但是用Delaunay插值拟合时，在边界处会产生一些奇异点；而这与之后进行的形态学操作有一些矛盾，因此，要对这些奇异点进行如下处理：用离奇异点相对近的点的灰度值来代替奇异点的灰度值，这一做法符合图像的基本特性且处理之后的效果显著。

2 试验结果

采用Matlab［7］编程，笔者提出的改进算法运行的时间为1．797s，直接用传统的经验模式分解运行的时间是6．797s，改进算法运行比传统方法快了5s；速度提高了73．56%，降低了时间复杂度，提高了运算效率。图像分解效果如图2所示。

3 结 语

对Norden Huang提出的EMD算法进行了改进，改进算法大大地降低了时间复杂度和空间复杂度。试验结果表明，该算法能大大的减小对计算机存储容量的要求，能有效的进行分解处理，在图像分解中有良好的效果，进一步研究可以应用在图像去噪、图像合成、图像复原等方面。

图2 试验结果

［1］Norden E H．A new method for nonlinear and nonstationary time series analysis：empirical mode decomposition and hillbert spectral analysis［J］．Proc SPIE，2000，4056：197-205．

［2］徐晓刚，徐冠雷，王晓通，等 ．经验模式分解 （EMD）及其应用 ［J］．电子学报，2009，03：151-153．

［3］冈萨雷斯 ．数字图象处理 ［M］．北京：电子工业出版社，2007：420-454．

［4］张合勇，任德明，赵卫疆，等 ．图像处理中二维经验模式分解的改进算法 ［J］．光学学报，2009，29（5）：1248-1253．

［5］Waston D F．Computing the n-Dimensional Delaunay Tessellation with Application to Voronoi Polytops［J］．The Computer Jour-nal，1981，24 （2）：16-20．

［6］Nunes J C，Guyot S，Delechelle E．Texture analysis based on local analysis of the bidimensional empirical mode decomposition ［J］．Machine Vision and Applications，2005，16 （3）：177-188．

［7］刘卫国．Matlab程序设计教程 ［M］．北京：中国水利水电出版社，2005．