基于形態學成分分析的合成孔徑雷達圖像去噪

王 燦 蘇衛民 顧 紅 邵 華

(南京理工大學電光學院，江蘇 南京 210094)

引 言

合成孔徑雷達(Synthetic Aperture Radar，SAR)雖然具有全天候全天時成像、高分辨和穿透性強等諸多優點，但是在SAR成像過程中，分辨單元的幅度受隨機相位的影響產生類似于噪聲的相干斑，相干斑嚴重影響了SAR圖像的質量[1]，使得人們對圖像目標的探測、分析和解譯能力降低.所以，相干斑抑制是SAR圖像處理的一個重要環節，引起了大家的廣泛關注.

在過去的十幾年里，很多相干斑抑制的方法被提出.傳統的方法分為兩類：多視處理和空域濾波處理.多視處理方法簡單，但卻犧牲了空間的分辨率.濾波處理方法根據像素周圍局部統計特性進行自適應空域濾波，能取得比較好的相干斑抑制效果，但是處理結果嚴重依賴濾波器的類型、方向和大小，因此很難較好的保留紋理和邊緣信息.

近幾年，隨著多分辨分析濾波技術的發展，基于小波變換理論的閾值去噪方法在SAR圖像相干斑抑制中得到應用[2-4].這些抑制方法相對傳統空域方法，能較好的實現區域平滑和紋理保留的折中，但是通用閾值去噪方法沒有考慮小波系數之間的相關性，所有小波系數統一被閾值處理，有過扼殺小波系數的傾向，圖像細節不能很好的保留.M.S.Crouse等人提出的基于空間域像素統計模型的去噪方法[5-8]，可以有效抑制SAR相干斑，但是對于大面積圖像，基于空間域模型的去噪方法受到限制.

2005年，Donoho提出的形態學成分分析方法MCA[9]通過不同的基對圖像稀疏表示可以把圖像的分段光滑部分和紋理部分分離出來.MCA方法主要應用于圖像修復和圖像分離[10-11].在文獻[9]中，MCA首次被用于相干斑抑制，但是因為它只考慮到圖像的平滑部分和紋理部分，在抑制噪聲的同時圖像的邊緣沒有得到很好的保留.本文基于一個新的超完備字典利用MCA思想進行相干斑抑制.首先，把圖像分為分段光滑部分，紋理部分，邊緣部分；然后，根據MCA思想，用一個包含脊小波(curvelet)的超完備的字典對圖像進行稀疏表示，把圖像的分段光滑部分、紋理部分、邊緣部分分離出來；最后利用分離出來的部分恢復圖像，從而可以有效抑制噪聲.同時因為圖像的特征信息被稀疏表示分離出來，這種方法可以很好的保留圖像的特征信息.實驗結果證明了本文算法的有效性.此外，本文給出了迭代算法中迭代步長c的上界，保證了迭代算法的收斂性.

1 SAR相干斑特性分析

SAR圖像相干斑被認為是一種乘積噪聲,有

Y(k，l)=F(k，l)X(k，l)，

(1)

式中：Y∈RN×N表示觀測圖像強度；X∈RN×N表示無噪聲圖像強度；F∈RN×N表示相干斑噪聲強度，k，l表示第k行，l列.對于L視的強度數據，F的概率密度函數服從均值為1方差為1/L的伽馬分布，即

(2)

式中Γ(·)表示伽馬函數.對式(1)兩邊進行取對數變換有

(3)

(4)

(5)

(6)

式中ψ(k，t)=(d/dt)k+1lnΓ(t)是第k階polygamma函數.通過對數變換及式(2)～(6)可知，乘性相干斑噪聲抑制問題可以轉化為加性高斯噪聲抑制問題.

2 基于MCA的SAR相干斑抑制

2.1 MCA

(7)

(8)

(9)

(10)

式中Φp，Φt，Φe依次記為可以稀疏表示平滑部分，紋理部分和邊緣部分的基.在本文中，選擇Daubechies3小波表示平滑部分的基Φp，離散余弦表示紋理部分的基Φt，curvelet小波表示邊緣部分的基Φe，變換域的系數為αp，αt，αe，并且這些系數是稀疏的.所以，利用稀疏表示把圖像的平滑部分，紋理部分和邊緣部分分離就稱為(MCA)[12]

(11)

2.2 迭代收縮算法SAR相干斑抑制

式(11)的求解是一個復雜的過程，所幸近幾年來，已經出現了很多算法，常用的有連續子空間最優化算法、快速迭代軟閾值算法和共軛梯度法[12-13].本文使用的是快速迭代收縮算法中的可分離代理函數法[14]，這種算法是由Daubechies 2004年提出的.為了表示方便，設

Φ=[Φp，Φt，Φe]；

(12)

α=[αp，αt，αe].

(13)

暫時忽略TV項.式(11)表示為

(14)

根據Daubechies的理論，添加下面一項為

(15)

式中，參數c的選擇必須保證函數dist是嚴格凸函數，也就是需要它的Hessian矩陣是正定的，即

cI-ΦTΦ>0，

(16)

式中，I表示單位矩陣.c必須滿足條件

(17)

式中λmax(·)表示最大特征值.結合式(14)和(15)，新的目標函數表示為

(18)

(19)

(20)

則式(19)再次簡化為

(21)

至此，式(18)的最優化問題最終簡化為酉形式的最優化問題[15].酉形式的最優化問題求解分為兩步，首先計算

(22)

利用收縮因子限制υ0的元素數目從而獲得估計值，其中收縮因子為

(23)

收縮因子把比較小的值映射到零，其他的值也向零收縮.最終得到α的最優估計值為

αopt=Sλ/c(υ0)

(24)

所以α的迭代公式為

(25)

把TV項考慮進去，分解的迭代公式如下：

(26)

式中：H是非抽樣Haar小波字典;Sγ表示以γ為閾值的收縮因子.與收縮去噪算法類似，利用冗余Haar小波收縮因子估計TV修正項.

根據上述討論，基于MCA的SAR圖像相干斑抑制算法的主要處理步驟如下：

(27)

設置初始圖像冗余部分為

(28)

(29)

3) 更新圖像冗余：

(30)

4) 轉到第(2)步，直到滿足收斂準則.

5) 最終輸出迭代估計稀疏系數為

(31)

6) 得到去噪后的圖像：

(32)

2.3 字典的選擇

選擇合適的字典是本算法的一個關鍵.字典必須能稀疏表示圖像的平滑部分、紋理部分和邊緣部分，它依賴于圖像的特征.根據多分辨分析的思想，小波變換的過程是圖像的分解過程，把圖像分了低頻部分和高頻部分，而非零系數主要集中在圖像的低頻部分，也就是圖像的平滑部分，所以選擇Daubechies3小波變換稀疏表示圖像的平滑部分.

離散余弦變換(Discrete Cosine Transform,DCT)的基是余弦函數，余弦函數有周期振蕩特性，因為紋理也具有周期變化的特點，所以本文利用DCT對圖像的紋理部分稀疏表示.

Donoho于1999年提出了Curvelet變換并且構造了Curvelet的緊框架，對于具有光滑奇異性曲線的目標函數，Curvelet提供了穩定的、高效的和近乎最優的表示.Curvelet變換的最大特點是各向異性，也就是Curvelet是一種具有方向性的基原子，可以對圖像多尺度多方向濾波，因此具有更強的表達圖像中沿邊緣信息的能力，所以本文選擇Curvelet稀疏表示圖像的邊緣部分.

2.4 參數c的邊界討論

(33)

即

(34)

在迭代算法中，固定的迭代步長引起迭代精度和迭代速度之間的矛盾.針對這個矛盾，本文選擇遞增幾何序列{ck}，它能同時提高算法性能和增加算法收斂速度.

2.5 閾值的選擇

相干斑抑制過程中，閾值的選擇尤為重要，它的選擇直接影響到去噪效果.如果閾值太小，閾值處理后的小波系數中包含過多的噪聲分量；如果閾值太大，將會丟失信號的一部分信息.針對這個矛盾，本文首先給出一個大的初始值λ/c，然后在每次迭代時逐次減少λ/c，直到減少到最小值，即

(35)

(a) 原始圖像 (b) Lee濾波后的圖像 (c) 小波閾值濾波

(d) 基于MCA平滑+紋理部分 (e) 本文方法 圖1 真實SAR圖像的濾波效果比較

3 原始數據的相干斑抑制實驗

在本節，為了驗證本文提出相干斑抑制算法的有效性，我們對實測機載條帶SAR圖像進行處理.實驗數據采用的是機載雷達X波段SAR成像的圖像數據，圖像是單視，大小為1 024×1 024像素，圖像的分辨率為1 m×1 m，圖像區域包括草地、農田、公路.分別利用Lee濾波，小波閾值去噪和基于MCA平滑加紋理算法和本文算法對SAR圖像進行去噪處理.其中Lee濾波算法中用到的鄰域窗口大小為5×5像素，小波域閾值去噪的閾值選擇為局部適應閾值.

圖1是SAR圖像經過幾種去噪算法后的效果圖.圖1(a)是原始SAR圖像，圖1(b)～(f)分別給出Lee、小波域軟閾值法、MCA方法小波加DCT和本文算法濾波后的圖像.可以看出Lee濾波器對噪聲起到一定的抑制作用，但是帶來圖像的分辨率降低.小波域軟閾值濾波效果好于Lee，但因為小波域軟閾值法對紋理、邊緣像素與噪聲像素統一處理，所以不可避免地破壞了一些紋理和邊緣結構.文獻[7]中的方法較小波域軟閾值法保留了圖像的紋理信息，但是部分邊緣信息丟失.本文方法利用curvelet小波稀疏表示SAR圖像的邊緣部分，提高了圖像變換域的稀疏度，因此不僅起到了相干斑抑制的作用，而且很好的保留了圖像中的紋理、邊緣等有用信息.如圖1方框中的局部放大圖圖2所示，因為原始SAR圖像中相干斑的存在，圖2(a)中微弱的周期信號不容易觀測到，利用已有算法對SAR圖像相干斑抑制后的圖像如圖2(b)～(d)所示，因為圖像的細節丟失，周期信號也丟失，而經過本文算法處理過的SAR圖像在抑制相干斑的同時，保留了圖像的細節信息，微弱的周期信號也可以清楚的觀測到，如圖2(e)所示.

(a) 原始圖像 (b) Lee濾波后的圖像 (c) 小波閾值濾波

(d) 基于MCA平滑+紋理部分 (e) 本文方法 圖2 真實SAR圖像的濾波效果的放大圖

表1列出通過不同算法濾波后圖像的等效視數和比值圖像的均值，由表1可知，本文方法在等效視數上是最高的，說明本文方法能有效地抑制同質區域的相干斑，并且由于我們使用超完備字典對SAR圖像的有用信息進行了稀疏表示，所以在圖像保持程度上也是最好的.

表1 不同算法的等效視數和比值圖像均值

4 結 論

利用curvelet能夠對圖像邊緣信息稀疏表示的特點，本文提出了一種基于形態學成分分析的SAR圖像相干斑抑制方法.該方法根據形態學成分分析把圖像分成平滑部分、紋理部分和邊緣部分，然后利用一個超完備字典中的不同基，通過迭代算法對圖像的平滑部分、紋理部分和邊緣部分分別進行稀疏表示，從而保留了圖像的有用信息，抑制了相干斑.實驗證明：與小波域等SAR圖像相干斑抑制方法比較，本文提出的算法在抑制相干斑的同時，提高了圖像邊緣等細節特征的保持能力，是一種穩定，有效的處理方法.另外，本文還給出了迭代算法步長的上界，保證了算法的收斂性.

[1] GOODMAN J W.Some fundamental properties of speckle[J].JOSA.1976，66(11):1145-1150．

[2] 陳 曦，張 紅，王 超.基于AOS非線性擴散的SAR圖像去噪研究[J].電波科學學報.2004，19(004):405-408．

CHEN Xi，ZHANG Hong，WANG Chao.A study of SAR images denoising based on AOS nonlinear diffusion[J].Chinese Journal of Radio Science，2004，19(004):405-408.(in Chinese)

[3] GAGNON L，JOUAN A.Speckle filtering of SAR images:A comparative study between complex-wavelet based and standard filters[C]∥Processing of The International Society for Optical Engineering. San Diego，1997．

[4] DONOHO D L.De-noising by soft-thresholding[J].IEEE Transactions on Information Theory，1995，41(3):613-627．

[5] VEGA M，MATEOS J，MOLINA R，et al.Bayesian TV denoising of SAR images[C]∥IEEE International Conference on Image Processing. Brussels，2011：165-168．

[6] PARRILLI S，PODERICO M，ANGELINO C V，et al.A Nonlocal SAR Image Denoising Algorithm Based on LLMMSE Wavelet Shrinkage[J].IEEE Transactions on Geoscience and Remote Sensing，2012，50(2):606-616．

[7] 吳 艷，王 霞，廖桂生.基于小波域隱馬爾可夫混合模型的 SAR 圖像降斑算法[J].電波科學學報.2007，22(2):244-250．

WU Yan，WANG Xia，LIAO Guisheng.SAR images despeckling based on wavelet and hidden Markov mixture model[J].Chinese Journal of Radio Science，2007，22(2):244-250.(in Chinese)

[8] CROUSE M S，NOWAK R D，BARANIUK R G.Wavelet-based statistical signal processing using hidden Markov models[J].IEEE Transactions on Signal Processing，1998，46(4):886-902．

[9] ELAD M，STARCK J L，QUERRE P，et al.Simultaneous cartoon and texture image inpainting using morphological component analysis(MCA)[J].Applied and Computational Harmonic Analysis，2005，19(3):340-358．

[10] BOBIN J，STARCK J L，FADILI J M，et al.Morphological component analysis:An adaptive thresholding strategy[J].IEEE Transactions on Image Processing，2007，16(11):2675-2681．

[11] STARCK J L，ELAD M，DONOHO D L.Image decomposition via the combination of sparse representations and a variational approach[J].IEEE Transactions on Image Processing，2005，14(10):1570-1582．

[12] BECK A，TEBOULLE M.A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J].SIAM Journal on Imaging Sciences.2009，2(1):183-202．

[13] PATEL V M，EASLEY G R，CHELLAPPA R.Component-based restoration of speckled images[C]∥ IEEE International Conference on Image Processing.Brussels，2011：2797-2800．

[14] DAUBECHIES I，DEFRISE M，DE MOL C.An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J].Communications on pure and applied mathematics，2004，57(11):1413-1457．

[15] ZIBULEVSKY M，ELAD M.L1-L2 optimization in signal and image processing[J].IEEE Signal Processing Magazine，2010，27(3):76-88．

[16] 徐 豐，金亞秋.多方位高分辨率SAR的三維目標自動重建(二)多方位重建[J].電波科學學報，2008,23(1):23-33．

XU Feng，JIN Yaqiu.Automatic reconstruction of 3D objects from multi-aspect part Ⅱ:Multi-aspect reconstruction[J].Chinese Journal of Radio Science.2008,23(1):23-33.