結合二維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分解與自適應高斯濾波的組合改進算法，將影像的高頻和低頻信息分開處理，對影像的高頻部分進行自適應高斯濾波去噪。通過兩個試驗的比較與評價，發現本文算法具有較大的峰值信噪比、平均梯度、結構相似性和較小的均方根誤差，且邊緣檢測結果也表明噪聲在被濾掉的同時，本文算法去噪后的影像能較多且更好地保留原始影像的細節信息和邊緣輪廓信息，取得較好的去噪效果，這對后期影像的使用具有重要的現實意義。