经验模式分解

经验模式分解

摘要

近些年来,随着计算机技术的高速发展与信号处理技术的不断提高,人们对图像的分析结构的要求也越来越高。目前图像处理已经发展出很多分支,包括图像分割、边缘检测、纹理分析、图像压缩等。经验模式分解(EMD)是希尔伯特-黄变换(Hilbert-HuangTransform)中的一部分,它是一种新的信号处理方法,并且在非线性、非平稳信号处理中取得了重大进步,表现出了强大的优势与独特的分析特点。该方法主要是将复杂的非平稳信号分解成若干不同尺度的单分量平稳信号与一个趋势残余项,所以具有自适应性、平稳化、局部性等优点。鉴于EMD方法在各领域的成功应用以及进一步的发展,国内外很多学者开始将其扩展到了二维信号分析领域中,并且也取得的一定的进展。但是由于二维信号不同于一种信号,限于信号的复杂性和二维数据的一些处理方法的有限性,二维经验模式分解(BEMD)在信号分析和处理精度上还存在一些问题,这也是本文要研究和改善的重点。

关键词：图像处理；信号分解；BEMD

Abstract

In recent years, with the rapid development of computer technology and the continuous improvement of signal processing technology, the demand for the analysis structure of the image is becoming more and more high. At present, many branches have been developed in image processing, including image segmentation, edge detection, texture analysis, image compression and so on. Empirical mode decomposition (EMD) is a part of Hilbert Huang transform (Hilbert-HuangTransform). It is a new signal processing method, and has made significant progress in nonlinear and non-stationary signal processing, showing strong advantages and unique analysis points. This method mainly decomposes the complex non-stationary signals into several single scale stationary signals with different scales and a trend residual term, so it has the advantages of adaptability, stationarity and locality. In view of the successful application and further development of EMD method in many fields, many scholars at home and abroad have expanded it to the two-dimensional signal analysis field, and have made some progress. However, because two dimensional signal is different from one signal, it is limited to the complexity of signal and the processing methods of two-dimensional data. Two-dimensional empirical mode decomposition (BEMD) still has some problems in the accuracy of signal analysis and processing, which is also the important point of research and improvement in this paper.

Key words: image processing; signal decomposition; BEMD

目录

摘要 (1)

第一章概况 (4)

2.EMD方法原理 (5)

2.1 本征模函数 (5)

2.2 .EMD分解过程 (5)

2.3.分解举例： (6)

3. BEMD分解原理 (8)

3.1 图像极值点的选取： (8)

3.2 Delaunay 三角剖分 (9)

3.3 基于三角网络的曲面插值 (11)

3.4 分解方法 (11)

3.5 BEMD 分解停止准则 (12)

4 二维经验模态分解在图像处理中的应用.................................. 错误！未定义书签。

4.1图像分解实例............................................................. 错误！未定义书签。

4.2图像降噪 ................................................................... 错误！未定义书签。5总结.................................................................................... 错误！未定义书签。参考文献................................................................................ 错误！未定义书签。

第一章概况

随着计算机技术的不断发展和其应用领域的不断扩展，数字图像处理技术得到了迅猛的发展，涉及信息科学、计算机科学、数学、物理学以及生物学等学科，因此数理及相关的边缘学科对图像处理科学的发展有越来越大的影响。近年来，数字图像处理技术日趋成熟，它广泛应用于空间探测、遥感、生物医学、人工智能以及工业检测等许多领域，并促使这些学科产生了新的发展。数字图像处理就是利用计算机对图像进行处理，其任务是将原图像的灰度分布作某种变换，使图像中的某部分信息更加突出，以使其适应于某种特殊的需求。目前数字图像处理已经发展出很多分支，包括图像分割、边缘检测、纹理分析、机器视觉等领域。多分辨率多尺度是人类视觉高效、准确工作的重要特征之一。自然产生的图像大多包含大量不同尺度的信息，这些信息在一幅图像中同时出现。而对图像的应用研究往往仅限于某一尺度或某些尺度上的现象，或者只需要某些尺度的信息；其他尺度的信息往往会对处理结果有不良影响，或者增大了处理的难度和复杂性。所以把图像信息按尺度进行分离非常必要。多尺度图像分解可以消除其他无用尺度信息对处理结果的影响，也简化了处理的难度和复杂性；也是图像目标识别和边缘检测等处理过程的预处理方法之一。

图像可以看作二维的随机信号，因此也可以通过对一维随机信号处理方法的二维扩展实现。经典的信号处理方法是傅立叶谱分析法，其必须对信号进行全局分析，不具备局部性，只适合线形和平稳信号分析。为了分析非平稳与非线性信号，人们发展了信号的时频分析法，如短时傅立叶变换、Wigner-Ville 分布以及小波变换等。由于所有的这些时频分析方法都是以 Fourier 变换为最终理论依据，故有一定难以克服的局限性，难以在时频分析的方法中取得突破。为了克服这些弊端，1998 年, Huang[1]等人提出了一种具有自适应时频分辨能力的信号分析方法——经验模式分解（EMD，Empirical Mode Decomposition）。鉴于经验模式分解在一维信号处理方面已经获得巨大的成功，将其推广到二维情况，将图像按尺度从小到大进行分离，小尺度信息包含了图像的细节信息，剩余的大尺度信息表达了图像的基本趋势和结构，将会给图像处理领域提供一种新的有效的数据处理手段。

2.EMD方法原理

大多数的工程实际信号都是由多频分量组成的,这种模式下的信号性质很复杂不容易分析。而单分量的信号在任一时间只有一个频率或一个频率窄带的信号,对人们来说分析起来即直观又明了,同时很容易进行进一步的深入处理。所以将多分量信号分解成若干个单分量的信号形式就是信号处理领域需解决的重点

EMD的分解过程就是将复杂信号单分量化,也就是说,对原始数据序列通过EMD分解成

若干个本征模函数(IMF)和一个残余分量。

2.1 本征模函数

在物理上,一个函数的瞬时频率有意义,那么这个函数序列必须是对称的、局部均值为零、具有相同的过零点和极值点数目。在这样的条件下,N.E.Huang等人提出了本征模函数(Intrinsic Mode Function,简称IMF)的概念[2]。假设任何信号都是由若干本征模函数组成,一个信号都可以包含若干个本征模函数,如果本征模函数之间相互重叠,便形成复合信号。

EMD分解的目的就是为了获取本征模函数。本征模函数形成条件有以下两个:

1.函数在整个时间范围内,局部极值点和过零点的数目必须相等或最多相差一个;

2.在任意时刻点,局部最大值的包络(上包络线)和局部最小值的包络(下包络线)平均值必须为零。

2.2 .EMD分解过程

EMD分解方法进行是需要一定条件

(1)整个信号中至少存在一个极大值点和一个极小值点。

(2)两个相邻的极大值点或者是两个相邻的极小值点间的时间间隔作为EMD分解的尺度特征。

(3)如果信号中不存在任何极值点,但是有拐点,可以通过一阶或者多阶微分求得极值点,最后再对分量求积分获得结果。

如果信号满足以上3个条件,则这个信号就可以进行EMD的分解。

具体的实现过程是：

设原数据序列为f(t)，首先对f(t)选取局部极值(包括极小值和极大值)，所有的极大值的插值曲线连接作为包络的上界，所有的极小值的插值曲线连接作为包络的下界，包络

的上下界将覆盖所有的数据，得到极大值包络数据E max(t)和极小值包络数据E min(t)。由这两条包络得到一个均值函数e1,即：

e1=E max t+E min(t)

2

e1与原始信号 f(t)的差量定义为分量d1，即

d1=f t?e1

考察d1是否为 IMF 函数。实际上d1局部包络均值不为零，也并未达到分离骑行波和使波形更加对称的目的。所以处理过程必须重复进行多次，直到得到内禀模式函数为止。在第二次的筛分过程中，d1被作为数据，

d1?e11=d11

其中e11为d1的平均值。

处理过程重复 k 次后，此时有：

d1(k?1)?e1k=d1k

定义c1=d1k，c1满足 IMF 条件，即为从原数据f(t)分离出的第一个 IMF。

SD=

|d1k?1t?d1k t|2

1k

2

T

t=0

通常 SD 值可以设置为 0.2～0.3，即满足 0.2

2.3.分解举例：

原始信号

极大值点

极大值包络曲线

极小值点

极小值包络包络

极大值和极小值的包络均值

从原始信号中减去均值后的余项

根据IMF满足的条件，考察其是否为IMF。若不是，则继续进行上面过程，直至得到IMF为止。

3. BEMD分解原理

BEMD的分解过程与一维的EMD分解类似,只是在极值点求取、包络面的构造等算法上存在不同,其次是二维的处理信息量很大,其空间和时间上的相关性也会影响到整个分解过程。信号可以进行BEMD分解的假设条件是[3]:

(1)在二维数据平面中,至少包含一个极大值点和一个极小值点;如果整个二维数据平面内不存在极值点,但可通过求导运算后求得一个极大值点和一个极小值点。

(2)二维分量的提取是以二维数据中的极值点之间的距离作为特征尺度。

3.1 图像极值点的选取：

图像的数据存储形式是栅格结构，栅格结构将图像划分成均匀分布的栅格，这些栅格就是通常所说的像素。图像数据数列中的局部极值点的求取常见的有两种方法: 一种方法是在像素点的8邻域中寻找极值点,即与8邻域点做灰度值比较确定极值点;

另一种方法是利用数学形态学方法寻找极值点[4]。

这里采用的是8邻域法寻找极值点。

对于一幅灰度图像,设左上角为坐标原点,水平方向为X轴,竖直方向为Y轴,建立二维坐标平面XOY,对于任意像素点(i,j)在图像中存在三种形式:(i,j)是图像的内部点、(i,j)是图像的边界点、(i,j)是图像的角点。如图所示。

图像内部点图像边点

图像角点

由于像素点的位置不同,对应周围邻域的大小不同,因此在提取灰度图像极值时,将像素点(i,j)的灰度值分别与对应邻域像素的灰度值进行比较。如果是图像内部的点, 将点

(i,j)的灰度值与周围的8邻域进行比较。同理,图像边界上的点和角点分别与周围5邻域

和3邻域进行比较,由于5邻域和3邻域是8邻域的特殊情况,因此这种提取灰度图像极值

的方法也称为8邻域比较法这种方法简单、检出的极值点完全，但计算量较大，需要时

间相对长。

3.2 Delaunay 三角剖分

找出图像数据序列中的所有极值点后,它们在二维平面上是一个散乱分布的状态,也就

是说它们构成了一个离散的点集,需要按一定的邻接关系将它们有序地组织起来,以便在空

间上进行曲面插值拟合,形成极值的包络面。因为三角形能够最好的逼近不规则的多边形曲面,因此常用三角形平面做逼近[5]。

Delaunay 三角剖分( DT, Delaunay Triangulation)可追溯到于 1907 年提出的 G. Voronoi 图，后来 Delaunay 在 1932 年提出了解决该剖分完整而实用的方法，近几十年

来该算法不断得到改进，以适应不同的应用。DT 也就是最近点意义下的 Voronoi图的直

线对偶图，Voronoi 图由许多胞元组成，每一个胞元包含点集中的一个点(每两点的垂直平分线将平面分成两部分，与某点对应的胞元实际是该点与点集中每个点垂直平分线在平面

内所形成的交集)，连接相邻胞元中两个点(直线段)便形成点集DT。

DT 具有很多性质，剖分中常用的性质有最小内角最大和最大空圆原则。

最小内角最大原则：对于一个凸四边形的两种剖分，DT 获得的两个三角形中的最小内角最大。

最大空圆原则：剖分中任一三角形的外接圆(三维为外接球面,高维为超球面)内不含有点集中的任何其他点。

最小内角最大原则又称为局部最优准则,这样获得的剖分称为局部最优(或局部等角)。

（1）（2）

两种不同的分割方式:

从图中可以看到，任何相邻四个孤立的点可以有两种不同的刨分形式，如果按照DT的最小内角最大的原则刨分的话，第一种形式满足条件要求，因为第二种形式的最小内角明显小于第一种形式的最小内角。第二种刨分模式会引起三角网格化出现大量的尖角形式的三角形式，所以要选择第一种刨分形式。

最大空圆原则

从图中可以看出这种剖分满足最大空圆原则,因为ABC的外接圆内不含有任何除了 ABC 三点外的其他离散点;同理,ACD也满足这样的要求。因为ABCD四点同时满足DT的两个条件,所以图中的这种剖分形式即为一个Delaunay三角剖分网格。将二维平面内所有的离散点都按照这两个条件进行DT三角剖分以后可得到整个平面内的三角网格化。

三角刨分基本步骤：

1、构造一个超级三角形，包含所有散点，放入三角形链表。

2、将点集中的散点依次插入，在三角形链表中找出外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，完成一个点在Delaunay三角形链表中的插入。

3、根据优化准则对局部新形成的三角形优化。将形成的三角形放入Delaunay三角形链表。

4、循环执行上述第2步，直到所有散点插入完毕。

3.3 基于三角网络的曲面插值

曲面插值在经验模式分解中占有重要地位,对于一维EMD分解来说,根据连续性和实时性等要求可选用不同的插值方法来满足要求,比如，当对分解精度要求不太高、对实时性要求高的情况下,可以采用线性插值或者最近插值等方法;如果在对分解精度要求高、对实时性要求不太高的情况下,则可以采用三次样条插值的方法,在保证所有数据的连续和平滑条件下,将所有数据特征都能有效的提取出来。对于二维经验模式分解来说,在选择插值方法方面同样值得着重考虑。

目前,在BEMD中较常见的插值方法有:径向基函数插值、线性插值、B-样条插值、单片样条插值等方法[6]。径向基函数插值需要求解大型的线性方程组,并且会产生大的矩阵数据,这样导致的结果是编程实现时会增大计算消耗以及存储空间需求很大,甚至会无法得到分解结果。线性插值虽然方法简单、运算量较小、计算较快,但是其插值效果过于线性,不能很好的模拟出理想的数据,会导致分解不完全的现象。B样条曲面插值的方法是使用B样条母线法,是使用“母线成面”的原理构造的。通过一个点集生成三次B样条曲线方程,然后得到n条B样条曲线,进而可得到曲面带,最后将曲面块连接到一起,形成一个插值曲面。

在这里用了三角域上的曲面插值,通过Ddaunay三角剖分得到整个平面上的三角网格图,然后对每个网格中使用多项式插值的方法可得到若干个三角曲面块,最后将这些三角曲面拼接到一起即可得到整个插值曲面,类似于分段函数中的各个函数段。

通过对极大值和极小值的插分曲面得到极大值曲面和极小值曲面。

3.4 分解方法

记f(x,y)为待处理的数字图像，其极大值包络曲面E max(x,y)和极小值包络曲面

E min(x,y)的代数均值记为E1(x,y)，即：

E1x,y=E max x,y?E min(x,y)

2

它与原图像数据f(x,y)的差值定义为D1(x,y)，即：

D1x,y=f x,y?E1x,y

D1(x,y)是f(x,y)的一个中间过程值，重复上述过程 k 次，直到D1k(x,y)是一个本征模函数，此时有：

D

1(k?1）

x,y?E1k x,y=D1k x,y

定义C1x,y=D1k x,y，则C1x,y即为分离出的第一个本征模函数。同样，必须确定一个准则用于每层筛分过程的停止。这可以由限制标准差的尺寸SD完成。判断第i层本征模函数筛分结束的标准差判别函数

一般取 SD ≤ 0.3。 然后把C 1 x,y 从原数据中分离得到余项R 1 x,y ，即

f x,y ?C 1 x,y =R 1 x,y

将R 1 x,y 作为新的数据，重复上述过程n 次，得到最终的表达式：

其中f x,y 为原图像数据，C i x,y 是分解后得到的较小尺度细节信息，R n x,y 是得到的最终大尺度趋势项。

3.5 BEMD 分解停止准则

二维经验模式分解的分解程序是一个二重循环过程，外层循环控制着整个图像信号分解的停止条件。在二维 EMD 分解过程中，如何判断分解过程是否应该结束呢？可以参考如下两个条件：

1大极小值点的总数小于 K （K 为一较小的正整数，可以人为设定）；

2 极大值点或极小值点都近似相等。

实际中，二维经验模式分解若满足两者之一，则强制分解过程结束。 若图像的 BEMD 分解停止条件太严格，分解得到的最后几个IMF 意义不大；若条件太宽松，可能丢失图像中的有用分量[9]. ∑∑==---=x x y y y x D k i y x D ik y x D k i SD 112

2

),()1()],(),()1([∑=+=n i n i y x y x y x f R C 1),(),(),(