结合二维EMD与自适应高斯滤波的遥感卫星影像去噪

王跃跃，陈 蓉,,于丽君，朱建峰，吴愈锋，陈炫炽

(1. 贵州大学矿业学院，贵州 贵阳 550025； 2. 贵州省非金属矿产资源综合利用重点实验室，贵州 贵阳 550025； 3. 中国科学院遥感与数字地球研究所，北京100101)

遥感卫星影像由于在接收和获取的过程中经常会不可避免地引入噪声，使得对后续遥感卫星影像的观察、研究和使用产生严重的影响。因此，在遥感影像数据的降噪处理中，既能保留原始影像细节信息，又能够去除影像中的噪声显得尤为重要。

传统的去噪方法主要分为两类[1]：一类是利用噪声和影像信息分布的频段不相同，通过算法将原始影像变换为不同频率的分量进行去噪，如文献[2]提出的双树复小波变换结合四阶偏微分影像去噪；另一类利用去噪滤波函数得到相应的模板，再用此模板对原始影像变换后灰度图像上的每一个像素的灰度值进行处理，如文献[3]采用张量多维滤波对影像进行去噪。但在实际处理中，前一种方法往往不能将影像的高频、低频信息分开处理，当去除噪声时会连同细节信息一起去除；后一种方法虽然在去除噪声的试验中效果还不错，但没有考虑用同一个模板去处理整个灰度图像中的不同像素的特性和适应性，特别是在处理纹理较多的影像时造成严重的模糊现象。近年来，随着EMD理论在信号和图像去噪方面研究的兴起，鉴于其在处理一维信号方面已经取得巨大的成功，将其推广到二维影像的处理并取得了一些成果，如文献[4—5]提出了一种二维EMD分解方法；文献[6]将经验模态分解运用于海洋遥感图像处理中，为其图谱分析和目标识别提供一种数据处理的新方法。其中也有一些改进算法，如文献[7]将图像转换为行向量、列向量和对角线向量后采用经验模式分解进行去噪；文献[8]针对图像利用随机微分对得到的多个固有模式函数和剩余函数进行滤波重组得到降噪影像。虽然这些方法在一定程度上取得了成功，但仍存在一些问题，如在处理噪声的同时未能将低频和高频分开，去噪时连同低频信息一起去掉，破坏了影像的基本趋势和框架，造成处理后的影像模糊等问题。因此本文在此基础上提出一种基于二维EMD与自适应高斯滤波的组合改进算法，通过两组试验表明，此方法可以用于遥感卫星影像的去噪，并且去噪后影像不仅保留了原始影像的大体趋势和基本框架，而且保留了原有的边界信息和纹理特征，有效地解决了使用传统去噪方法去噪后导致影像模糊的问题。

1 遥感影像的EMD去噪理论

对于一幅影像f(x,y)进行二维EMD分解，将其自适应分解为IMF分量，即固有模态函数，所得到的固有模态函数有以下特征：

(1) 无论什么影像数据都可以分解为IMF分量，并且任意两个及以上的IMF分量都可以合成为新的影像数据。

(2) 每一个固有模态函数的上下包络线关于时间轴局部对称，局部极大值和极小值与时间轴的交点数目相同，如不同最多相差为1。

基于以上特征，固有模态分解可以将其他算法不能处理的非线性、非平稳的影像数据自适应地分解为多个固有模态函数[9]，但是其分解过程必须满足以下条件：①影像数据平面中必须要有极值点，即至少一个极大值点和一个极小值点；②在时域的性质中极值距离具有重要的作用；③当影像数据中无极值点，仅包含拐点时，可对拐点的微分方程求解来获取极值点，并对求解结果进行一个类似“筛选”的过程[10-11]，即可得到极值点。二维EMD分解实现流程如下：

(1) 计算上下包络线的极值点和平均值。求取影像序列的所有极大值点U1(t)和极小值点U2(t)并对其进行拟合，从而算出一个关于上下包络线的平均值L1(t)，即

(1)

(2) 二维EMD进行影像分解。为了获取经验模式分解后的分量P1(t)，将原影像序列中每个点的值X(t)分别减去上下包络线的平均值L1(t)。

P1(t)=X(t)-L1(t)

(2)

当P1(t)满足前文提到的3个条件，那么分解出来的第一阶IMF即P1(t)。如不满足前文3个条件，则把P1(t)当作新的原始影像序列，并重复步骤(1)、(2)n次，直至P1n(t)=P1(n-1)(t)-L1n(t)满足条件，则P1n即新的第一阶IMF，把它记为z1(t)。

z1(t)=P1n(t)

(3)

(3) 二维EMD分解结束条件。定义原始影像序列为q(t)，令

F(t)=q(t)-z1(t)

(4)

式(4)得到的F(t)与原影像数据极为相似，对F(t)重复以上步骤，便可得到每一阶IMF分量，当在分解过程中出现的分解结果满足假设条件时结束二维EMD分解。

2 遥感影像自适应高斯滤波去噪

图像的去噪往往需要计算图像中每一个像素中心点附近(1+2k)×(1+2k)大小的领域Sxy[12]，这样即可用其来构造自适应高斯滤波器。而传统的高斯滤波器一般均表示为[13-14]

(5)

式中，σ为标准差；k为高斯核都是预先设定的。这样的高斯滤波器往往在对图像进行去噪时会连同原始影像信息一起去掉，从而造成影像模糊。对于自适应高斯滤波器，标准差和高斯核都应该是作为式(5)的函数，即

(6)

其中最关键的一步就是建立Sxy和σxy、kxy之间的对应关系，通过前人的总结和查阅文献得出三者的关系为

Sxy=σxy=kxy

(7)

最终将式(6)求出的结果离散化为(1+2kxy)×(1+2kxy)大小的模板，即为自适应高斯滤波的处理模板，用其去处理前文计算出来的IMF分量图即可得出自适应高斯滤波去噪的结果图。

3 遥感影像的二维EMD与自适应高斯滤波去噪

任何一幅影像都是由不同的频率成分组成的，影像的高频部分代表着原始影像数据的大量细节信息和噪声，低频部分代表着影像的趋势信息[15-19]。因此，本文在此基础上，先将影像进行低通滤波处理，得到高频和低频两部分，低频部分保留不变，对高频信息进行二维EMD分解，再把每个分量进行自适应高斯滤波去噪重建原始图像，这样既保留了高频部分的细节信息，又保留了原始影像的低频信息，并可以防止原始影像的趋势信息丢失。本文基于Matlab R2014a软件平台将影像高低频分开，对高频部分采用二维EMD分解和自适应高斯滤波相结合的算法，对传统的EMD算法和高斯滤波算法进行改进。具体改进算法的基本流程如图1所示，以下介绍本算法的详细步骤。

(1) 获取高频影像。在Matlab R2014a中编程将原始遥感影像进行空间域滤波(低通滤波)获取低频部分和高频部分。由于原始影像通过空间域滤波后灰度值变小，导致原始影像的高频部分对比度和细节信息不明显，因此对其进行线性拉伸，拉伸的尺度为原始影像的灰度范围。

(8)

式中，H输入为输入影像；G低通滤波(*)为空间域滤波中的低通滤波；K低频为低频影像；K高频为高频影像；X线性拉伸为线性拉伸；T拉伸后为线性拉伸后的高频影像。

(2) 计算拉伸前后高频影像的均值

(9)

式中，求出的a1为未进行线性拉伸前的均值；a2为进行线性拉伸后的均值。

(3) 二维EMD分解。高频影像进行拉伸后作经验模式分解，得到4个不同频率的固有模态函数，即

Vi=P(K高频)i=1,2,3,4

(10)

式中，Vi为IMF分量；P(*)为EMD分解。

(4) 对每个IMF分量进行高斯自适应去噪。由于EMD分解后每个IMF分量图的灰度值和像素值都不一定相同，因此采用自适应高斯滤波对每一个分量图进行处理可以更好地保留每个分量图的结构特征与细节信息，下文有详细概述，这里不再赘述。

Zi=F(Vi)i=1,2,3,4

(11)

式中，Zi为自适应高斯滤波处理后的IMF分量；F(*)为高斯自适应去噪。

(5) 自适应高斯滤波处理后的IMF分量重构高频影像

(12)

式中，X为自适应高斯滤波处理后的IMF分量重构的高频影像。

(6) 获取结果影像。将自适应高斯滤波处理后的IMF分量重构的高频影像和原始影像空间域滤波后得到的低频影像进行叠加获取结果影像。其中由于前文对原始影像的高频部分进行过线性拉伸，为了让结果影像在处理前后的灰度值保持不变，必须计算灰度的均值差，具体处理方式见式(4)—式(6)。

J=K低频+X-(a2-a1)

(13)

式中，J为结果影像；K低频为原始影像通过空间域滤波后得到的低频影像。

4 试验结果与分析

4.1 试验1

为了验证和对比本文改进算法的有效性，在试验1中将二维EMD分别结合维纳滤波、高斯滤波、中值滤波和小波进行去噪。具体算法流程如图2所示。

截取贵州省花溪区一景SPOT卫星5 m高分辨率全色波段影像的局部区域影像为例，大小为502×502像素，灰度级为256，现对该清晰影像加入方差为0.1的高斯噪声。在试验中使用二维EMD方法将图像分解为3个IMF分量及一个剩余函数图像，加入噪声的影像图和分解后的分量图如图3所示。

本文试验中将加噪影像通过试验1算法和本文算法，并作了效果对比。试验1算法与本文算法的区别在于对图3得到的IMF分量采取处理的滤波函数不同，试验结果见图4所示。

从图4各种方法去噪影像的目视效果来看，试验1中所用算法和本文算法都能较好地去除噪声，但试验1中所用算法与本文算法去噪效果比较而言还有一些不足，如二维EMD结合维纳滤波去噪后的图4(a)中，其去噪效果并不均匀，影像的大部分噪声得以去除，但局部还是存有较多噪声；二维EMD结合高斯滤波去噪后的图4(b)整个影像还残留少量噪声，影像的轮廓和细节信息也不够清晰；二维EMD结合小波去噪后的图4(d)去噪效果最差，不仅影像还存有一些噪声，轮廓也不清晰；二维EMD结合中值滤波去噪后的图4(c)与图4(a)、(b)和(d)的去噪效果相比有所提升，但在图4(c)的中间部位和边缘部位仍存在轮廓不清晰、细节信息不突出、部分影像纹理模糊不见的问题；而经过本文算法的去噪影像图4(e)，不仅噪声去除明显，影像的边界、轮廓清晰可见，原始影像的细节纹理信息也得到较好保留。

4.2 试验2

为了验证和评价试验1中算法与本文算法的去噪效果，本文从以下两个方面来展开：

(1) 6个评价指标：由于试验1中有原始影像作为参考对比影像，因此在试验1中采用峰值信噪比(PSNR)、均方根误差(MSE)、均值、熵、平均梯度和结构相似性(SSIM)作为评价指标去评价去噪后的影像，见表1。

表1 试验1算法及本文算法去噪后结果影像质量评价指标统计表

由表1可以看出，本文算法得出的评价指标与试验1中各算法得出的评价指标都相差不大，说明去噪后的影像没有失真，这也可以由表中的结构相似性得以证明，本文算法的结构相似性相比试验1中的4种组合算法略高(结构相似性的最大值为1，最小值为0)，说明经过本文算法去噪后的影像与原始图像有更高的吻合度和相似度。而本文算法的熵只小于试验1算法之二维EMD结合高斯滤波去噪的熵，熵越大说明图像所含信息越丰富，此为本文考虑每一个IMF分量图中像素的特异性而采用高斯滤波来作IMF分量的自适应去噪的原因。本文算法的平均梯度相比试验1中的4种组合算法较大，说明遥感影像经过本文算法降噪后具有清晰的边界轮廓和细节信息。通过计算去噪后的影像和原始不含噪声影像间的峰值信噪比和均方根误差得出本文算法去噪后的峰值信噪比最大、均方根误差最小，这说明本文算法在去噪效果上优于试验1中的4种组合算法。

(2) 影像边缘检测：遥感影像去噪的目的是为了获得更清晰的影像，以便后续的处理与应用。遥感影像去噪后的效果好坏，一个关键性的评价因素为去噪后的影像是否保留了原始影像的边界轮廓信息。因此本小节将对试验1中各种组合算法及本文算法去噪后的影像进行边缘轮廓检测，边缘检测采用Matlab中自带的canny算子，图5为检测结果。

从图5中可以看出，在经过canny算子边缘检测后，本文算法去噪后的影像边缘轮廓检测结果最好，二维EMD结合高斯滤波其次，二维EMD结合中值滤波第三，二维EMD结合维纳滤波和二维EMD结合小波检测结果最差。从检测结果图的总体来看，本文算法得出的图5(e)线条比其他4种要密集，二维EMD结合维纳滤波和二维EMD结合小波得出的图线条最为稀疏，说明在去噪过程中原始影像的一些细节信息被当作噪声去掉，未能保留原始影像的边界轮廓信息；而二维EMD结合高斯滤波和二维EMD结合中值滤波去噪后影像的边缘检测图中局部边界轮廓信息丢失， 且线条不连续多呈分叉

状，说明在去噪过程中原始影像的结构受到破坏和去噪后影像有噪声存在影响了边缘检测的结果；本文算法检测结果图中线条密集、连续，特别是在原始影像最右边的模糊区，经过本文算法去噪后都能大致检测出完整的边界，由此说明本文算法去噪后残留的噪声少，原始影像的结构在去噪过程中被破坏程度小，保留了较多的细节信息。

从以上两个方面的阐述可以看出，经过本文算法去噪后的影像，无论是在前文的6个评价指标上还是在边缘检测后的效果上都明显优于试验1中的算法。

5 结 语

遥感卫星影像的去噪在影像处理领域中一直是一项关键技术。本文将二维EMD分解理论引入遥感影像的去噪中，不仅拓宽了其应用范围，而且对于成像复杂、易受外界自然因素干扰的遥感卫星影像，EMD去噪理论展示了它相对于传统算法能更好地处理非平稳、非线性信号的优势。本文提出的二维EMD分解与自适应高斯滤波的组合改进算法，将影像的高频和低频信息分开处理，对影像的高频部分进行自适应高斯滤波去噪。通过两个试验的比较与评价，发现本文算法具有较大的峰值信噪比、平均梯度、结构相似性和较小的均方根误差，且边缘检测结果也表明噪声在被滤掉的同时，本文算法去噪后的影像能较多且更好地保留原始影像的细节信息和边缘轮廓信息，取得较好的去噪效果，这对后期影像的使用具有重要的现实意义。