一种改进的BEMD阈值去噪算法在单目视觉定位中的应用

柳 笛，陈熙源，方文辉

（东南大学 仪器科学与工程学院 微惯性仪表与先进导航技术教育部重点实验室，南京 210096）

图像作为视觉定位系统解算位姿的一种重要载体。图像获取和传输过程中会受到外在或内在因素的影响，从而产生噪声污染。这些噪声降低了图像的质量，致使图像中的重要信息丢失，给通过图像解算载体位姿信息带来很大的困难。因此图像的降噪处理和稀疏表达，是视觉定位分析的重要基础[1]。

近些年来，一些基于分解和重构的算法被应用于图像去噪，例如各种滤波器、傅里叶变换和一些较复杂的超小波等算法（小波去噪[2]、剪切波变换去噪[3]、曲波变换去噪[4]）。这些方法根据噪声的高频特性，在对图像信号进行分解和重构时，通常将变换系数较大的高频信号作为噪声而滤除，从而达到图像去噪的目的。事实上，图像的纹理和边界等细节同样包含了大量的高频信息，因此图像中的细节信息往往会和噪声同时被平滑。小波去噪方法能够在多方向上进行图像去噪，然而它们需要预先设计分解基底，如果分解达不到最优则会造成不能对图像进行彻底降噪或者将一些有用的细节信息作为噪声而滤除。1990年，美国国家宇航局的Huang等人[5]提出一种用于非线性和非平稳数据分析的自适应时域分解方法-经验模态分解（EMD），该方法具有很强的局部自适应分解特性，在一定程度上克服了小波分解算法的不足并在一维信号分析中取得了广泛的应用。Park等人[6]采用基于分位数的经验模式分解（QEMMD）对二维图像进行了去噪研究，该方法通过将含有噪声的低阶IMF去除来达到降噪的目的，但其忽略了低阶本征模态函数中仍含有的少量有用信息，从而可能造成原始图像信号的变形。

基于以上去噪算法的不足，为能够将低阶本征模态函数中噪声去除的同时并保留其含有的少量有用信息，本文提出了一种改进的BEMD阈值去噪算法，并对仿真的含噪声图像和单目视觉定位系统获取的图像进行了降噪实验。结果表明改进后的BEMMD阈值去噪算法较小波去噪和BEMD-DTT去噪方法的降噪性能有明显改善，并显著提高了单目视觉定位系统的定位精度。

1 二维经验模态分解原理

经验模态分解（EMD）是基于信号局部特征时间尺度，将非线性、非平稳信号分解成不同尺度本征模态函数的一种方法。二维经验模态分解（BEMD）正是由一维 EMMD拓展而来，并可用于图像处理分析。通过BEMD分解图像得到的本征模态函数（IMF）依然要满足两个条件：1）在整个信号数据集上，极值点和过零点的数目必须相等或者至多相差一个；2）在任何一个时间点上，图像像素的局部极大值点和局部极小值点构成的上下包络面的平均值为零。对于一个二维图像I(x,y)进行BEMD分解，其具体步骤如下[7]：

1）初始化图像，令Iri,jj-1(x,y)=I(x,y)作为输入图像，且i=1,j=1 。

2）计算出Irri,j-1(x,y)的极大值点和极小值点，根据极值点构造出二维图像的上包络面emaxi,,j-1(x,y)和下包络面emini,j-1(xx,y)，并求取其平均曲面：

3）将输入图像减去平均包络面得到：

4）判断hi,j-1(x,y)是否满足 IMF的约束条件SDD＜ε。如果满足，则第i个本征模态函数IMFi(x,y)为hhi,j-1(x,y)；如果不满足，则j=j+1，令Iri,j-1(x,y)=hi,j--2(x,y)并返回步骤2）直至经过k次循环满足约束条件，从而得到IIMFi(x,y)。其中ε一般取值小于0.25，

5）求hi,j-1(x,y)取残余分量。i=i+1,Iri,0(x,y)=Iri-1,0(x,y)-IMFFi-1(x,y)，如果得到的残余分量中极值点的个数多于两个，则将Iri,,0(x,y)作为新的图像转至步骤2）继续完成分解；否则分解结束，得到不能再分解的残余分量Irkk(x,y)，于是原图像可以表示为：

本文采用以上步骤完成了对测试图像I(x,y)（BBaboon图像）的分解，其中测试图像加入了均值为0，方差为0.01的随机噪声，其分解结果如图1所示，图11 的(b)～(m)为测试图像的12个本征模态函数图像，(n)为测试图像的残余分量。

图1 加噪后图像 I( x,y)的BEMD分解：(a) 原始的噪声图像；(b)-((m) 第1至第122阶本征模态函数；(n) BEMD的残余项Fig.1 BBEMD of the noiisy image I( x, y);(a) original nooisy image;(b)-(mm) represent the first to twelfth IMMFs;(n) residual oof BEMD

2 改进的BEMD阈值去噪算法

BEMD是自适应的将图像分解成不同尺度特性的IMF。目前常用的BEMD图像降噪方法主要包括两种：一种是根据高频噪声主要集中在低阶本征模态函数中的特点，将前几层低阶模态函数舍弃而达到降噪的目的，这种方法虽然简单，但它忽略了低阶模态函数中含有的少量细节信息；另一种是结合其它去噪算法（如小波去噪）对每层本征模态函数去噪后重构，该方法能在一定程度上达到降噪效果，但是对不含噪声的本征模态函数进行降噪处理显然增加了去噪的耗时性，同时也是非必要的。针对以上常用方法存在的不足，本文提出了一种改进的BEMD阈值去噪算法。

2.1 噪声主导模态函数与信号主导模态函数的区分

图像在经过 BEMD分解后得到一组本征模态函数，该组模态分量可分为两类，一类是同时含有噪声和图像信息的噪声主导模态分量，另一类是只含有有用信息的信号主导模态分量。为了能够最大程度的抑制图像中的噪声，准确区分出噪声主导模态函数和信号主导模态函数显得尤为重要[10-11]。

本文提出了一种基于2ℓ准则和概率密度函数（PDF）的二维图像噪声主导模态函数与信号主导模态函数区分方法。该方法首先使用标准正态核密度函数对本征模态函数与原图像的概率密度函数进行估计；然后通过几何距离对每个本征模态函数的概率密度函数与原图像概率密度函数之间的相似性进行测量。根据贝叶斯对概率密度函数的解释，概率密度函数反映的不仅仅是信号的频率，它还体现了信号的状态信息，概率密度函数包含有关其对应IMF的完整信息，因此可以通过概率密度函数相似性度量来识别包含有信号主要特征的IMF。在模态区分的过程中，假定J作为噪声主导模态函数与信号主导模态函数的边界索引值。原始图像与本征模态分量的概率密度函数之间的相似性可通过下式求出：

其中，PDFS(⋅)表示相似性测量值，PDF(⋅)表示概率密度函数，dist(⋅)表示概率密度函数之间的几何距离。

本文中选用欧几里得距离来进行概率密度函数相似性的测量，其被定义为：

式中，P和Q表示概率密度函数。

在所计算的一组欧几里得距离中，其值逐渐增加，达到局部最大值后开始减小时所对应的模态分量就是第一个信号主导模态函数，因此噪声主导模态函数与信号主导模态函数的边界索引值J被定义为：

图2中展示的是图1中各本征模态分量的概率密度函数与测试图像的概率密度函数。

图2 测试图像I(x,y)与其各本征模态函数(IMF1～IMF12)概率密度函数Fig.2 PDFs of the imageI(x,y) and its corresponding IMF(IMF1～IMF12)

从图2中可以看出，从第一阶本征模态函数开始，其概率密度函数与测试图像概率密度函数之间的相似性逐渐降低，直到第9阶本征模态函数达到最低，其后两者的相似性突增并趋于平稳。

图3所示的测试图像与其各本征模态函数间的概率密度函数相似性测量值更能直接的反映出这种趋势关系，因此测试图像的噪声主导模态函数与信号主导模态函数的边界索引值J为10。

图3 测试图像I(x,y)与其各本征模态函数间的概率密度函数相似性测量值Fig.3 PDF-based similarity versus IMFs of the image(,)I(x,y)

基于BEMD分解的传统去噪方法，如BEMD-DT方法在区分噪声主导模态函数与信号主导模态函数时，根据噪声的高频特性，认为噪声主要集中在前两阶IMF中，因此这种方法认为前两阶IMF是噪声主导模态函数，剩余的IMF是信号主导模态函数。噪声主导模态函数与信号主导模态函数区分的精确与否将直接影响图像去噪的效果。为了证明本文所提区分方法的优越性，第3节将BEMD-DT算法与本文所提算法的去噪效果进行了对比。

2.2 噪声的压缩

由于在信号采集过程中引入的噪声等原因造成了EMD分解过程中存在位置敏感性误差问题，从而使分解得到的IMF丧失了具体的物理意义。为了解决该问题，WU[12]等人提出了一种总体平均经验模态分解法（Ensemble Empirical Model Decomposition,EEMD），该方法可以有效消除位置敏感性误差问题，但是EEMD存在运算量大，运算速度慢等缺点。

为了克服EEMD方法的缺点，崔冰波[13]等基于白噪声零均值的特点，提出了一种将白噪声加入含噪信号中使更多的高频噪声压缩到低阶IMF中的方法，该方法减少了低阶IMF的分解误差，显著提高了EMD的分解精度，并且该方法只对初步估计为噪声占优的IMF分量进行噪声辅助数据分析，所以在一定程度上减小了计算量，提高了运行速度。将以上几种应用于一维信号分解的改进方法应用到二维图像分解中，在消除位置敏感性误差方面可以取得较好的效果，但是仍存在计算量大、运算速度慢的缺点。

针对以上问题，本文提出的改进方法通过多次改变噪声分布位置重新构建含噪图像，对重构的多幅含噪图像进行求平均处理，并以较小的计算量将图像噪声压缩至低阶IMF中，以达到减小低阶IMF分解误差，显著提高BEMD分解精度的目的。

该改进方法的具体描述如下：

步骤①：对含噪图像I(x,y)进行传统的BEMD分解，得到L阶IMF；

步骤②：对I(x,y)与各阶IMF的概率密度函数相似性进行测量，得到噪声主导模态函数与信号主导模态函数的边界索引值J，并对前1J-个IMF进行重构，得到信号

步骤③：使用小波对信号IP(x,y)进行降噪处理，信号IP(x,y)降噪后表示为IP(x,y)；

步骤④：计算图像信号IP(x,y)中存在的实际噪声In(x,y)，In(x,y)=IP(x,y)-Ip(x,y)，并完成对最后L-J+1个IMF、残余分量Ir(x,y)及IP中所含有用信息的重构：

步骤⑤：任意改变噪声序列In(x,y)的位置分布，得到新的噪声分布序列Ia(x,y)：

步骤⑥：构建含有不同噪声分布并与原图像具有相同信噪比的图像信号：

步骤⑦：重复执行步骤⑤⑥K次，得到K幅不同噪声分布的图像I1(x,y),I2(x,y),…,IK(x,y)；

步骤⑧：将K幅不同噪声分布的图像信号进行平均处理，得到噪声压缩后的图像信号：

步骤⑨：使用BEMD方法对得到的噪声压缩图像信号I(x,y进行分解，得到N阶IMF。

在上述步骤⑦中，关于重复迭代次数K的设置，从理论上来说，K的值越大，含噪图像信号的分解越准确，但此时会降低图像分解的速度。在此处选取K为20进行噪声压缩实验，仅为了证明噪声压缩方法的有效性。关于迭代次数K取不同值时对图像去噪效果的影响将在第3节进行统计分析。

根据上述步骤对含噪图像进行噪声压缩，图4中展示了测试图像I(x,y)在噪声压缩前后的对比结果，可以看出噪声压缩后的图像质量较原图像略有改善。将噪声压缩图像进行 BEMD分解并对本征模态函数与噪声压缩图像的概率密度函数进行估计。

图4 噪声压缩前后图像对比：(a) 原噪声图像 I( x, y))；(b) 声压缩图像(x,y)Fig.4 Image before and after noise ccompression:(a) original nooisy image I(x,y);((b) noise comprr essed image ( x, y)

图5 噪声压缩图像(xx, y)与其各本征模态函数(IIMF1～IMF12)的概率密度函数Fig.5 PDFs of the noise-commpressed image( x, y)and its corrresponding IMMF(IMF1～IMF12)

图6 噪声压缩图像( x, y)与其各本征模态函数间的概率密度函数相似性测量值Fig.6 PDFF-based similarrity versus IMFFs of the noisse-compressed iimage ( x, y)

为了与EEMD算法进行对比，我们假定EEMDD 算法的总体平均次数和本文所提方法的重复迭代次数均取值为K(K=1,2,…)。因此，一次完整的EEMD过程需要进行K次BBEMD分解和一次IMF的总体平均处理。本文所提噪声压缩方法包括两次BEMMD分解和一次对不同噪声分布的图像信号进行平均处理的过程。

图7中展示了在不同K值下，本文方法和EEMDD方法压缩噪声所消耗的时间。通常，EEMMD需要几十甚至上百次的总体平均次数才能取得较好的噪声压缩效果。因此，在这种情况下，本文算法的执行速度明显优于EEMD算法。

图7 本文方法和EEMD方法在不同KK值下压缩噪声所消耗时间的对比Fig.7 Commparison of thee time consumedd by the proposed metthod and EEMDD method in compressed noise at differeennt K values

图8和图9分别展示了在K取值为40情况下，EEMD分解I(x,y)图像得到的各阶IMF与I(x,y)之间的概率密度函数相似性测量值和各阶IMF的概率密度函数分布。从图9可以看出，EEMD分解结果较本文算法增加了一阶IMF。从图8可以看出，与图6相比，EEMD算法将噪声压缩到前7阶IMF中，较本文压缩算法提高了一阶。因此，在噪声压缩速度和效果上，本文算法均优于EEMD算法。

图8 原噪声图像I(x,y)与通过EEMD得到的各阶IMF之间的概率密度函数相似性测量值Fig.8 PDF-based similarity versus IMFs of the original noise imageI(x,y)obtained by EEMD

图9 EEMD方法分解得到的各阶IMF与原噪声图像I(x,y)的概率密度函数Fig.9 PDF of original noise imageI(x,y) and IMFs obtained by EEMD

2.3 阈值去噪

根据上述噪声压缩图像分解结果，本文将使用软阈值法[14]对噪声主导模态函数进行降噪处理。该方法首先估计出各阶噪声主导模态函数的一个阈值，如果噪声主导模态函数中某点的幅值大于该阶函数的阈值，则保留该点，否则认为该点为噪声并将该点幅值置0。含噪图像I(x,y)经改进 BEMD分解重构后得到噪声压缩图像(x,y)，并且(x,y)被分解为N阶IMF和一个残余项Ir(x,y)，将(x,y)与各阶IMF的概率密度函数进行相似性测量得到噪声主导模态函数与信号主导模态函数的边界索引值J，对(x,y)进行阈值降噪处理得到图像IQ(x,y)。降噪后图像IQ(x,y)可表示为：

式中，Ti表示IMFi中的阈值。

3 仿真与应用分析

3.1 仿真分析

为了验证本文算法的有效性和适应度，对加入不同程度随机噪声的 Baboon图像进行去噪实验，并与小波阈值降噪、BEMD-DT降噪方法进行对比，其中小波阈值降噪选择了两种小波基（Sym4与Db6）。

为了能够比较客观的评价本文所提算法的去噪效果，比较同类降噪方法的优劣，本文中使用了均方误差（Mean Squared Error,MSE）和峰值信噪比（Peak Siignal to Noise Ratio,PSNR）两项指标对各类算法的降噪效果进行评价，均方误差和峰值信噪比的表达式如下：

1）均方误差：

式中，m×n表示图像大小，I′(x,y)是原无噪声图像，IQ(x,y)为降噪后的图像，均方误差越小表明去噪效果越好。

2）峰值信噪比：

式中，R表示图像的灰度量级，其取值为 255，峰值信噪比越高说明经过去噪处理后的图像失真越小。

在对含噪图像进行去噪前，为了选定合适的迭代次数K，使用本文算法对含有均值为零、方差为0.01随机噪声的Baboon图像在K取不同值时的去噪效果进行了统计分析。从图10中可以看出，随着迭代次数的增加，图像的峰值信噪比也相应提高。图11则表明随着迭代次数的增加，图像去噪所消耗的时间也随之增大。在综合考虑图像分解速度和分解精度的情况下，在以下的实验中K均取值为40。

图10 迭代次数与所对应图像去噪后峰值信噪比Fig.10 AAchieved PSNRof denoised image with respect to the number of iterations

图11 迭代次数与图像去噪所消耗时间Fig.11 Timespent on imagee denoising withh respect to the number off iterations

在实验中，分别对含均值为0，方差为 0.005、0.01、0.05和0.1随机噪声的Baboon图像采用小波阈值、BEMD-DT以及本文方法进行了去噪实验。图12中展示了对含有标准方差σ为0.01随机噪声的Baboon图像的去噪实验结果。从图12中可以看出，虽然小波和BEMMD-DT去噪算法的去噪效果达到了令人可以接受的程度。但与本文算法相比，其去噪效果还存在一定差距。为了定量分析这种差距的大小，表1给出了含有随机噪声的图像去噪前后的PSNR和MMSE。

从表中可以看出在不同小波基下的图像去噪效果是不同的，因此在使用小波阈值降噪方法时，小波基以及小波分解层数的选取至关重要，它们将影响最终图像的去噪效果。在以上几种降噪算法的去噪结果中，传统的BEMMD-DT算法的去噪效果最差，其原因在于该算法在图像去噪过程中并没有准确区分出噪声主导模态分量和信号主导模态分量，而是只对前两阶IMF进行了阈值去噪。本文算法对含噪图像去噪后获取的PSNR均高于其它两种算法，并且随着噪声强度的增加峰值信噪比提高的越明显。

图12 随机噪声的Baboon图像去噪结果对比：(a) Baboon图像 ；(b) 加噪图像；(c) sy4方法；(d) Db6方法；(e)BEMMD-DT方法；(f) 本文方法Fig.12 Comparison of denoising results of the different methods for the Baboon image with random noise:(a) original Baboon image;(b) noisy imagge;(c) Sym4;(dd) Db6;(e)BEMD-DT;(ff)) proposed

表1 不同去噪方法对含有不同方差随机噪声Baboon图像去噪前后的MSE和PPSNR（dB）Tab.1 MSE and PSNR before and after denoiising of Baboon image containing different random noise(dB)

3.2 在单目视觉定位中的应用分析

在室外实验过程中，我们使用如图13所示的单目视觉定位系统（其主要传感器参数见表2和表3）获取小车在直线和转弯（非直线运动）两种运动情况下定位所需的图像以及IIMU数据，并分别对使用本文算法、小波算法和BEMDD-DT算法降噪处理后的图像进行定位解算，其中Symm4和Db6小波分解层数均为3。在定位过程中，采用文献[15]中所提算法构建相似变换模型对小车轨迹进行解算，并将通过IMU数据解算的轨迹作为基准。

图13 单目视觉定位系统Fig.13 Monocularvisual positioning system

表2 IMU设备主要性能参数Tab.2 Main performance parameters of IMU

表3 工业相机主要性能参数Tab.3 Main performance parameters of industrial camera

图14为小车在室外环境下进行直线运动和转弯（非直线运动）运动后解算出的运动轨迹。从图14中可以看出，通过本文算法处理后解算出的小车位置最接近于小车的真实轨迹，小波次之，BEMD-DT方法取得的定位效果最差。进一步对比分析图15-17可知，通过改进BEMD阈值去噪算法处理并解算出的北向位置误差、东向位置误差和航向角误差均小于其它几种方法得到的误差。造成这种现象的原因是本文算法不仅有效去除了图像中的噪声，还较完整的保留了定位所需要的图像细节信息。

图14 位置解算结果：(a) 直线运动情况下位置解算结果；(b) 非直线运动情况下位置解算结果Fig.14 Positioncalculation ressults:(a) Position calculation reesults under lineear motion;(b)Position calculation results under non-lineear motion

图15 东向位置误差：(a)直线运动情况下东向位置误差；(bb) 非直线运动情况下东向位置误差Fig.15 East posittion error:(a)East position errror under linear motion;(b) East position error under non-linnear motion

图16 北向位置误差：(a) 直线运动情况下北向位置误差；(b) 非直线运动情况下北向位置误差Fig.16 North position error:(a) North position error under linear motion;(b) North position error under non-linear motion

图17 航向角误差：(a) 直线运动情况下航向角误差；(b) 非直线运动情况下航向角误差Fig.17 Heading angle error:(a) Heading angle error under linear motion;(b) Heading angle error under non-linear motion

表4 直线运动情况下不同方法解算结果的RMSETab.4 RMSEs of different methods in linear motion

为了定量分析本文算法、小波降噪算法（Sym4和Db6）和BEMD-DT降噪算法获得的状态量误差，表4和表5中分别统计了小车在直线和非直线运动情况下的东向位置、北向位置和航向角的均方根误差（RMSE）。从表中可以看出，本文算法显著提高了单目视觉定位系统的鲁棒性。在直线和非直线运动情况下，本文算法较Sym4小波去噪算法关于东向位置误差、北向位置误差和航向角误差分别平均改善了74%、64%、54%；较Db6算法分别平均改善了60%、48%、39%；较BEMD-DT算法分别改善了73%、96%、84%。

表5 非直线运动情况下不同方法解算结果的RMSETab.5 RMSEs of different methods in non-linear motion

4 结 论

针对单目视觉定位系统获取的图像含有较大噪声的问题，本文中提出了一种改进BEMD阈值去噪算法并将其用于解决该问题。该算法能准确区分出噪声主导模态函数和信号主导模态函数，在减小噪声对低阶IMF分解精度影响的同时可实现将噪声压缩到低阶 IMF中，其在有效去除噪声主导模态函数所含噪声的过程中还可以很好地保留其含有的少量图像细节信息。用此方法对含随机噪声的Baboon标准图像进行处理，取得了较好的降噪效果。将本文算法应用于单目视觉定位系统，并与小波去噪算法（Sym4和Db6）和BEMD-DT去噪算法相比，在东向位置、北向位置和航向角误差方面，本文方法较其他方法均有大幅度改善。