一種改進的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去噪算法相比，在東向位置、北向位置和航向角誤差方面，本文方法較其他方法均有大幅度改善。