利用經(jīng)驗?zāi)B(tài)分解和主成分分析的SAR圖像相干斑抑制

王文波，趙 攀，張曉東

1.武漢科技大學(xué) 理學(xué)院，湖北 武漢 430065；2.北京師范大學(xué) 遙感科學(xué)國家重點實驗室，北京 100101；3.武漢大學(xué)測繪遙感信息工程國家重點實驗室，湖北 武漢 430079；4.衛(wèi)星海洋環(huán)境動力學(xué)國家重點實驗室，浙江 杭州 310012

1 引 言

合成孔徑雷達（synthetic aperture radar，SAR）圖像會受到其自身固有的相干斑噪聲的影響，相干斑噪聲嚴(yán)重影響圖像的質(zhì)量，給圖像的自動解譯造成很大的困難［1］。因此，相干斑抑制成為SAR圖像處理的關(guān)鍵步驟，也是后繼SAR圖像特征提取、分割、識別的基礎(chǔ)。為了有效去除相干斑噪聲，國內(nèi)外學(xué)者進行了大量的研究，分別提出了基于核獨立成分分析［1］、基于線性多通道最優(yōu)濾波求和［2］以及基于小波分解［3－4］等多種 SAR圖像濾波算法。利用經(jīng)驗?zāi)B(tài)分解（empirical mode decomposition，EMD）進行相干斑抑制為SAR 圖像去噪提供 了 一種新的思路［5－6］，EMD 是文獻［7］提出的一種新的時頻分析工具，此方法將信號分解成有限個內(nèi)蘊模態(tài)函數(shù)（intrinsic mode function，IMF）和一個均值趨勢函數(shù)之和。EMD通過篩選的迭代程序來取得基底IMF，IMF會隨信號的不同而改變，它是一組完全數(shù)據(jù)驅(qū)動的自適應(yīng)基底，不需要事先設(shè)定基底和分解層數(shù)。近年來，EMD算法被逐漸應(yīng)用到了SAR圖像的相干斑抑制中，并取得了很好的去噪效果。文獻［8—10］對含乘性噪聲的SAR圖像進行4個方向的EMD分解，根據(jù)邊緣所在的方向從IMF中選擇相應(yīng)數(shù)據(jù)進行部分重構(gòu)，有效地去除了相干斑噪聲；文獻［11—12］應(yīng)用EMD分解對相干條紋噪聲進行了有效去除；文獻［13］利用EMD對SAR的干涉相位圖像進行分解，然后通過部分重構(gòu)的方法有效去除了干涉相位圖像中的相干斑噪聲；文 獻 ［14］應(yīng) 用 復(fù) 數(shù) EMD（complex EMD，CEMD）同時分解SAR圖像實部和虛部，根據(jù)閾值確定噪聲項IMF，對非噪聲項IMF進行部分重構(gòu)以去除窄帶相干噪聲，也取得了很好的去噪效果。

在文獻［8—14］的EMD相干斑抑制算法中，都是采用部分重構(gòu)的方法進行。在部分重構(gòu)去噪算法中，對EMD分解后的前幾項IMF認為全部由噪聲構(gòu)成直接去除，對剩下的IMF進行累加重構(gòu)，以達到去噪的目的。文獻［15］研究表明，含噪信號經(jīng)EMD分解后，前幾項IMF中仍含有較多的細節(jié)信息，而且在后面的IMF中也分布有不同程度的噪聲，因此采用直接部分重構(gòu)的方法對SAR圖像去噪時，會損失大量的細節(jié)信息，噪聲也不能得到徹底的去除。為了克服EMD部分重構(gòu)去噪算法的不足，本文采用PCA對SAR圖像的IMF作進一步處理，提出了一種改進的EMD部分重構(gòu)去噪算法。首先利用PCA對SAR圖像的IMF進行分解，然后根據(jù)IMF中所含噪聲的能量自適應(yīng)選擇合適的主分量重構(gòu)IMF，進一步去除IMF中的噪聲；最后對處理后的所有IMF進行累加重構(gòu)，得到去噪SAR圖像。

2 經(jīng)驗?zāi)B(tài)分解

3 基于PCA的EMD部分重構(gòu)相干斑抑制算法

3.1 imfk中所含噪聲的能量計算

將SAR圖像進行對數(shù)變換并進行均值調(diào)整后，乘性相干斑噪聲近似變?yōu)榱憔档募有愿咚拱自肼暎?6］。設(shè) 對 數(shù) SAR 圖 像 為I0，設(shè)I0經(jīng)EMD分解后，第k層IMF為imfk＝［d1d2…dN］，定義imfk的能量ε（imfk）為

式中，N表示imfk的長度。為了方便討論，令fk＝imfk，假設(shè)fk＝gk＋Vk，其中g(shù)k表示fk中所包含的信號成分；Vk表示fk中的噪聲成分。設(shè)Vk的能量為ε（Vk），根據(jù)含噪信號經(jīng)EMD分解后的噪聲能量模型［17］，可對imfk中噪聲的能量ε（Vk）進行近似計算

式中，β≈0.719；ρ≈2.01。但簡單地令ε（V1）＝ε（imf1），不考慮imf1中所含的細節(jié)信息，會導(dǎo)致噪聲能量估計偏大，不利于信號細節(jié)信息的保留。本文首先利用PCA對imf1進行去噪，再估計各層IMF中噪聲的能量。

3.2 利用PCA對各層IMF進行去噪

對數(shù)SAR圖像經(jīng)EMD分解后，各層IMF中都含有不同程度的噪聲，希望利用PCA進一步去除各層IMF中的噪聲。記

則CX＝E｛（X－mX）（X－mX）T｝為fk的協(xié)方差矩陣。設(shè)λ1≥λ2≥…≥λN是CX的特征值，φ1、φ2、…、φN是對應(yīng)的特征向量，記Φ＝［φ1φ2…φN］T，則Φ是一個正交矩陣。定義Y＝［Y1Y2…YN］T＝Φ（X－mX），根據(jù)PCA的分解特性可知，X中的噪聲分布在所有的成分分量Yi中，而信號主要集中在前幾層的成分分量中。如果選擇前H個特征向量構(gòu)成新的變換矩陣ΦH＝［φ1φ2…φH0 … 0］T，定義新的變換YH＝ΦH（X－mX），則通過反變換，可以由YH得到原信號X的近似值~X

3.2.1 利用PCA對imf1去噪時主分量的選擇方法

由 EMD對噪聲信號的分解特性可知［17－18］，imf1基本由噪聲構(gòu)成，只含有少量的信號細節(jié)信息，該類信號經(jīng)PCA分解后，信號基本只集中在第一個成分分量中［19－20］。因此本文在利用PCA對f1進行去噪時，僅保留第一個主成分分量進行重構(gòu)。即令X＝f1，X利用PCA去噪后的信號為＝［φ10 0 … 0］Y＝φ1Y1，此時從X中所刪除的噪聲為

從而可求出f1中所含噪聲的能量ε（V1）＝ε（ΔX）。根據(jù)EMD分解中噪聲能量分布模型式（2）和ε（V1）的值，可對fk（k≥2）中所含噪聲的能量進行計算。

3.2.2 利用PCA對imfk（k≥2）去噪時主分量的選擇方法

以下討論利用PCA對fk（k≥2）進行去噪時，主成分分量的選擇方法。如果在PCA重構(gòu)時，通過選擇合適的H，使得式（3）中刪除的噪聲ΔX的能量與fk中所包含噪聲Vk的能量相同，亦即在選擇H時，使得ε（ΔX）＝ε（Vk），則可認為噪聲被全部去除，保留下來的主成分分量就是理想的不含噪聲的信號gk。為了方便計算并減小誤差，對上式進行一定的變形，上式等價于

由式（3）可知，所刪除的噪聲ΔX的能量為

而信號X的能量為

3.3 基于PCA的EMD相干斑抑制的步驟

本文提出的聯(lián)合PCA和EMD進行相干斑抑制的步驟如下。

步驟1 計算SAR圖像I的邊緣方向［8］。

步驟2 對I進行對數(shù)變換和均值處理，設(shè)處理后的對數(shù)SAR圖像為IL。

步驟3 對IL按照4個方向分別進行EMD分解，對各分解結(jié)果按以下方式進行去噪：① 令X＝fk，求X協(xié)方差矩陣CX的特征根λ1≥λ2≥…≥λm和對應(yīng)特征向量｛φi｝；② 對于f1，取H＝1，進行重構(gòu)去噪得到，根據(jù)式（4）計算f1中噪聲的能量ε（V1），并根據(jù)公式（2）估計fk（k≥2）中的噪聲能量ε（Vk）；③ 對于fk（k≥2），根據(jù)式（8）、式（15）確定H值，進行重構(gòu)去噪得到；④ 令，則即為imfk去噪后的結(jié)果，對（k≥1）進行累加重構(gòu)，得到去噪后的對數(shù)圖像。

步驟4 對4個方向去噪后的對數(shù)圖像進行指數(shù)變換，結(jié)合邊緣方向計算最終去噪圖像Id。

4 試驗結(jié)果及分析

為了客觀評價SAR圖像的相干斑抑制效果，本文采用以下3個標(biāo)準(zhǔn)對去噪結(jié)果進行評價［21］：① 等效視數(shù)（ENL）；② 均值比MR；③ 邊緣保持指數(shù)ESI。為了驗證本文算法的有效性，試驗中采用兩幅SAR圖像作為測試圖像：英國Bedfordshire地區(qū)，X波段，2視SAR圖像，大小為256像素×256像素，見圖1（a）；來自美國Sandia國家實驗室的3視X－SAR圖像，大小為512像素×512像素，見圖2（a）。在試驗中，分別采用傳統(tǒng)的EMD部分重構(gòu)去噪算法 （conventional EMD partial reconstruction，EMD－CPR ）［17］、改 進 的EMD閾值去噪算法 （clear iterative EMD interval－thresholding，EMD－CIIT）［15］和本文方法（improve EMD partial reconstruction，EMDICPR）進行相干斑抑制。在EMD－ICPR算法中，對各層IMF進行PCA處理時，采用以下方式：對IMF進行分塊PCA分解；采用文獻［19］中的方法進行訓(xùn)練樣本集選擇，試驗中去噪模塊的尺寸設(shè)為6×6，樣本集訓(xùn)練模塊的尺寸設(shè)為30×30。圖1（a）和圖2（a）由不同算法去噪后的結(jié)果分別如圖2（b）～圖2（d）和3（b）～圖3（d）所示。

圖1 SAR圖像1相干斑抑制結(jié)果Fig.1 Speckle reduction comparison of SAR image 1

圖2 SAR圖像2相干斑抑制結(jié)果Fig.2 Speckle reduction comparison of SAR image 2

從圖1和圖2的視覺效果上分析可知，EMD－CPR方法去噪后，邊緣保持較好，比較清晰，但有大量噪聲未能去除而且有較多的細節(jié)信息損失；EMD－CIIT方法的去噪效果要好于傳統(tǒng)部分重構(gòu)算法，在噪聲去除和細節(jié)保持兩方面都有較好改進。EMD－ICPR方法得到的結(jié)果整體清晰度更高，與EMD－CIIT的去噪結(jié)果相比可以看出：EMD－ICPR方法在平滑區(qū)域噪聲去除得更徹底；在邊緣和紋理細節(jié)區(qū)域細節(jié)損失更少，較多的細小邊緣都能夠被保留，而EMD－CIIT的去噪結(jié)果中，一些細小的邊緣和紋理無法得到有效保持。為了定量對比各個算法的抑制平滑效果，在圖1（a）和圖2（a）中分別選擇3個同質(zhì)區(qū)域，圖1（a）中所選區(qū)域的大小為25像素×25像素、35像素×35像素、50像素×50像素，圖2（a）中所選區(qū)域的大小為59像素×93像素、46像素×72像素、25像素×33像素，計算去噪后這些區(qū)域的ENL，試驗結(jié)果由表1所示。由表1的試驗數(shù)據(jù)可知，EMD－ICPR方法對同質(zhì)區(qū)域相干斑的抑制能力最優(yōu)，與EMD－CPR算法相比ENL平均提高約60.474，與EMD－CIIT相比平均提高約22.462。

試驗結(jié)果的均值比MR和邊緣保持指數(shù)ESI分別由表2和表3給出。從表2和表3中的數(shù)據(jù)可知，EMD－ICPR方法去噪后，圖像的均值比也是最優(yōu)，很好地保持了圖像的輻射特性，與EMDCPR和EMD－CIIT相比，MR分別平均降低了0.136和0.076；而且EMD－ICPR方法去噪后的SAR圖像，在水平和垂直方向的邊緣保持指數(shù)ESI均高于另外兩種方法，與EMD－CPR相比，ESI平均提高了0.122 8，與EMD－CIIT相比，邊緣保持指數(shù)ESI平均提高0.056 7。這說明本文方法在消除噪聲的同時，有效地保持了圖像的邊緣細節(jié)信息。通過以上分析可以看出，不論是從視覺效果，還是從定量的客觀評價指標(biāo)來看，EMD－ICPR方法都具有一定的優(yōu)勢，與現(xiàn)有EMD去噪算法相比，有較好的改進。

表1 均勻區(qū)域等效視數(shù)比較Tab.1 ENL comparison of homogeneous regions

表2 SAR圖像1均值比和ESI的比較Tab.2 MR and ESI comparison of SAR image 1

表3 SAR圖像2均值比和ESI的比較Tab.3 MR and ESI comparison of SAR image 2

5 結(jié) 論

EMD分解是一種完全數(shù)據(jù)驅(qū)動的自適應(yīng)多尺度分解算法，與傳統(tǒng)的多尺度分解如金字塔或小波分解不同，EMD不需要預(yù)先給定基地和分解層數(shù)，而是根據(jù)數(shù)據(jù)本身特性自適應(yīng)地確定基地和分解層數(shù)，特別適用于非線性、非穩(wěn)態(tài)數(shù)據(jù)的分析與處理。本文分析和研究了對數(shù)SAR圖像經(jīng)EMD分解后各層IMF中噪聲的分布特性和噪聲能量的估計方法，給出了利用PCA對各層IMF進行去噪處理時主成分分量的選擇方法，在此基礎(chǔ)上提出了一種改進的EMD相干斑抑制方法。本文方法在一定程度上克服了EMD部分重構(gòu)去噪算法細節(jié)信息損失較大且噪聲不能完全去除的缺點，通過兩組SAR圖像進行了試驗比較，結(jié)果表明本文提出的利用PCA和EMD進行相干斑抑制的方法具有良好的去噪效果。由于EMD分解是一個迭代計算的過程，計算復(fù)雜度較高，因此本文提出的相干斑抑制算法的計算效率較低，在今后的研究中如何采用并行算法提高計算效率是一個重要問題。此外，由于EMD方法本身的理論還不成熟，如何從理論上證明對數(shù)SAR圖像的IMF中噪聲仍為加性噪聲，并結(jié)合試驗構(gòu)造IMF系數(shù)比較精確的分布模型，也是今后EMD相干斑抑制中需要重點研究的問題。

［1］ ZHANG Zhongshan，YU Jie，YAN Qin，et al.Research on Polarimetric SAR Image Speckle Reduction Using Kernel Independent Component Analysis ［J］.Acta Geodaetica et Cargographica Sinica，2011，40（3）：290－295.（張中山，余潔，燕琴，等.基于核獨立成分分析的極化SAR圖像相干斑抑制［J］.測繪學(xué)報，2011，49（3）：290－295.）

［2］ YANG Shenbing，LI Bingbai，SHEN Shuanghe，et al.Structure Retaining Linear Multi－channel SAR Image Speckle Filter［J］.Acta Geodaetica et Cargographica Sinica，2006，35（4）：364－370.（楊沈斌，李秉柏，申雙和，等.基于特征保持的線性多通道最優(yōu)求和SAR圖像濾波算法.測繪學(xué)報，2006，35（4）：364－370.）

［3］ ZHANG Jun，LIU Jian.A Speckle Reduction Algorithm by Soft－thresholding Based on Wavelet Filters for SAR Images［J］.Acta Geodaetica et Cargographica Sinica，1998，27（2）：119－124.（張俊，柳健.SAR圖像斑點噪聲的小波軟門限濾除算法［J］.測繪學(xué)報，1998，27（2）：119－124.）

［4］ STIAN S，ELTIFT T.A Stationary Wavelet－Domain Wiener Filter for Correlated Speckle［J］.IEEE Transactions on Geosciences and Remote Sensing，2008，46（4）：1219－1230.

［5］ ELGAMEL S A，SORAGHAN J.Empirical Mode Decomposition－Based Monopulse Processor for Enhanced Radar Tracking in the Presence of High－Power Interference［J］.LET Radar，Sonar and Navigation，2011，5（7）：769－779.

［6］ OLUFEMI A，VLADIMIR A，AUROOP R.Empirical Mode Decomposition Technique with Conditional Mutual Information for Denoising Operational Sensor Data［J］.IEEE Sensors Journal，2011，11（10）：2565－2575.

［7］ HUANG NE，SHEN Z，LONG SR，et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis［C］∥Proceedings of the Royal Society of London.London：［s.n.］，1998：903－995.

［8］ HAN C M，GUO H D，CHANG L W，et al.A Novel Method to Reduce Speckle in SAR Images［J］.International Journal of Remote Sensing，2002，23：5095－5101.

［9］ HAN Chunming，GUO Huadong，WANG Changlin.Speckle Suppression Using the Empirical Mode Decomposition［J］.Journal of Remote Sensing，2002，6（4）：266－271.（韓春明，郭華東，王長林.利用經(jīng)驗?zāi)B(tài)分解方法抑制SAR斑點噪聲［J］.遙感學(xué)報，2002，6（4）：266－271.）

［10］ HAN Chunming.The Study of SAR Image Speckle Filter［D］.Beijing：Graduate University of Chinese Academy of Sciences，2003.（韓春明.SAR 圖像斑點濾波研究［D］.北京：中國科學(xué)院研究生院，2003.）

［11］ MESSINA A，VITTAL V，HEYDT G.Nonstationary Approaches to Trend Identification and Denoising of Measured Power System Oscillations［J］.IEEE Transactions on Power Systems，2009，20（4）：1798－1807.

［12］ BERNINI B，GALIZZI E，F(xiàn)EDERICO A，et al.Evaluation of the 1DEmpirical Mode Decomposition Method to Smooth Digital Speckle Pattern Interferometry Fringes［J］.Optics and Lasers in Engineering，2007，45：723－729.

［13］ YUE H Y，GUO H D，HAN C M，et al.A SAR Interferogram Filter Based on The Empirical Mode Decomposition Method［J］.Geoscience Remote Sensing Symposium，2001，5：2061－2063.

［14］ FENG Zhou，XING Mengdao，BAI Xueru，et al.Narrow－Band Interference Suppression for SAR Based on Complex Empirical Mode Decomposition［J］.IEEE Geoscience and Remote Sensing Letters，2009，6（3）：423－427.

［15］ KOPSINIS Y，MCLAUGHLI S.Development of EMDBased Denoising Methods Inspired by Wavelet Thresholding［J］.IEEE Transactions on Signal Processing，2009，57（4）：1351－1362.

［16］ HUA X，LELAND E P，F(xiàn)AWWZA T U.Statistical Properties of Logarithmically Transformed Speckle［J］.IEEE Transactions on Geoscience and Remote Sensing，2002，40（3）：721－727.

［17］ FLANDRIN P，RILLING G，GONCALVES P.Empirical Mode Decomposition as a Filter Bank［J］.IEEE Signal Processing Letters，2004，11（2）：112－114.

［18］ WU Z，HUANG N E.A Study of the Characteristics of White Noise Using the Empirical Mode Decomposition Method［C］∥Proceedings of the Royal Society of London.London：［s.n.］，2004：1597－1611.

［19］ ZHANG Lei，RASTISLAV L，WU Xiaolin.PCA－based Spatially Adaptive Denoising of CFA Images for Singlesensor Digital Cameras［J］.IEEE Transactions on Image Processing，2009，18（4）：797－812.

［20］ MURESAN D D，PARKS T W.Adaptive Principal Components and Image Denoising［C］∥Proceedings of International Conference on Image Processing.［S.l.］：IEEE，2003：14－17.

［21］ HUANG Shiqi，LIU Daizhi.Research on Method and Application of Speckle Noise Reduction of SAR Image［J］.Acta Geodaetica et Cartographica Sinica，2006，35（3）：245－250.（黃世奇，劉代志.SAR圖像斑點噪聲抑制方法與應(yīng)用研究［J］.測繪學(xué)報，2006，35（3）：245－250.