基于2D-PWVD的SAR圖像相干斑噪聲抑制

劉 聰，李言俊，劉保元

(1.中國人民解放軍駐西飛公司軍事代表室， 西安710089;2.西北工業大學航天學院， 西安710072)

0 引言

由于合成孔徑雷達(SAR)圖像的成像特點，相干斑噪聲是SAR圖像的固有特性，為實現對SAR圖像的準確解譯，研究SAR圖像的噪聲抑制非常有意義。目前，已提出的SAR圖像相干斑噪聲抑制方法主要分為時域濾波方法和變換域濾波方法兩類。時域濾波方法主要包括:Lee濾波算法、Kuan濾波算法、Frost濾波及其增強算法、Gamma-MAP濾波算法等。變換域濾波算法主要是基于小波變換的濾波算法。

基于時域濾波的所有方法都是利用SAR圖像的統計特性，利用一個滑動窗口與原圖像做卷積，達到對圖像平滑的目的。這類方法實際上是以犧牲圖像細節為代價，達到噪聲抑制的目的，這就造成了SAR圖像斑塊抑制效果與圖像的邊緣細節保留之間的矛盾，通常增大滑動窗口的尺寸，能提高斑點噪聲抑制的效果，但同時也造成更多圖像結構信息的丟失，為保證二者的兼顧，窗口尺寸的合理選擇非常困難。

基于小波變換的相干斑噪聲抑制方法是將斑點噪聲轉換為高頻各層次的小波系數，根據高頻子域系數特征設計不同的處理策略，達到相干斑噪聲抑制的目的。相比于時域濾波方法，基于小波變換的濾波方法能比較好地實現區域平滑和紋理保留的折中。圍繞高頻子域系數特征的處理策略，出現了多種基于小波變換的噪聲抑制方法。Donodo、Johnstone［1］等提出了稱為“小波收縮”的信號去噪方法，在此基礎上，Devore與 Lucier［2］、B.Sanker［3］等提出了大量優化方法［4］。Romberg［5］基于多種隱馬爾科夫模型(HMM)，利用小波系數間的相關依存關系，對小波系數進行了重新估計，使其逼近原始圖像的小波系數，達到噪聲抑制的目的。

利用偽魏格納分布對圖像進行分解能得到圖像的低頻和高頻能量譜圖，一般認為圖像的邊緣噪聲能量主要集中在高頻區域，受小波軟閾值SAR圖像噪聲抑制方法的啟發，結合魏格納分布的特性，本文闡述了基于二維偽魏格納分布的SAR圖像噪聲抑制方法。

1 二維魏格納分布理論

Wigner于1932年首先提出了Wigner分布的概念，并應用于量子力學領域。1948年Ville將Wigner分布引入信號分析領域，應用于信號的檢測和信號細節分類。因此Wigner分布又稱為Wigner-Ville分布(WVD)，Wigner-Ville分布屬于雙線性時頻分析Cohen＇s類的特殊形式，具體表達式為

式中:*代表復數的共軛;t為時間分量;ω為頻率分量;W(t，ω)為信號f(t)的魏格納分布。

Wigner-Ville分布具有一系列好的性質，如奇偶性、虛實性、時間邊緣性和頻率邊緣性等。Jacobson和Wechsler在20世紀80年代將Wigner-Ville分布拓展到二維信號和三維信號的處理中，并闡述了多維Wigner-Ville分布與一維Wigner-Ville分布具有相同的特性［6-8］。二維 Wigner-Ville分布(2D-WVD)定義如下

式中:Wf(x，y，u，v)為信號 f(x，y)的二維維格納分布;(x，y)為時間分量;(u，v)為頻率分量。

離散形式定義為

信號的WVD變換存在交叉項的問題，交叉項是信號的N元成分經WVD變換后形成的N(N-1)/2項幅度，為各信號成分兩倍的干擾項，干擾項的存在嚴重影響對信號的準確認識，偽魏格納分布(PWVD)通過在時域加窗能有效抑制多元成分信號的交叉項，達到正確認識信號的目的，二維PWVD(2D-PWVD)定義為

由式(5)可以看出，選擇大小為(2N1-1)×(2N2-1)大小的時域窗口，針對二維信號的每個像素點進行2D-PWVD，可以得到整個二維信號的PWVD變換圖。

2 圖像的二維偽魏格納分解

圖像為二維數據，根據表達式(3)，對M×N圖像進行2D-WVD變換，需要建立一個M×N×M×N大小的四維數組，若再對M×N×M×N四維數組進行相關處理，計算量相當大，不利于現實應用，且圖像多元成分的交叉項影響也相當嚴重。根據局部頻譜思想，本文采用表達式(5)對原圖像進行2D-PWVD變換。在時域上，選取大小為N1×N2(要求N1=N2，且為奇數)的窗口函數在原圖像進行滑動，得到原圖像每個像素點大小為N1×N2的頻譜圖像。定義I(i，j)為數字圖像空間坐標(i，j)處的像素值，Xij為 I(i，j)的根據式(5)計算得到的2D-PWVD，其中

N1×N2為二維滑動窗口大小。代表像素I(i，j)在(k，j)頻率點的PWVD 值，其中k=1，2，…，N1，l=1，2，…，N2。提取圖像各空間位置的‖‖，按照圖像對應空間位置組成“能量譜圖”，一幅圖像能得到N1×N2幅與原圖像同樣大小的“能量譜圖”。分解流程如圖1所示。分解圖像的排列方式如圖2所示。

圖1 圖像數據的2D-PWVD分解流程

圖2 分解圖像的排列方式

選取時域窗口函數h(k，l)=1(矩形窗)，大小為5×5，通過圖1的變換重組流程，將lenna原始圖像分解為25個不同的頻段圖(圖3)。從分解重組圖可以看出，圖3(11)相對于lenna原始圖像，圖像邊緣模糊，可以理解為是原始圖像經過2D-PWVD變換重組的低通濾波圖像，分解圖像(12)～(55)是原始圖像邊緣的信息體現，可以理解為是原始圖像經過2D-PWVD變換重組的高通濾波圖像。仔細分析lenna圖像的分解圖像，可以看出:(12)和(15)、(13)和(14)、(21)和(51)、(31)和(41)、(22)和(55)、(23)和(54)、(24)和(53)、(25)和(52)、(32)和(45)、(33)和(44)、(34)和(43)、(35)和(42)分別是完全相同的。因此，lenna圖像經過2D-PWVD變換重組，實際上形成了13個分解圖像。推廣到所有時域窗口的情況下，分解重組后形成的圖像數為(N1×N2+1)/2。

圖3 lenna圖像的2D-PWVD 5×5分解圖

在圖3的高頻分解圖像中，(12)和(13)對原圖像的水平向變化劇烈的邊緣敏感，(21)和(31)對原圖像的垂直向變化劇烈的邊緣敏感，(22)和(33)對原圖像反對角線和對角線方向變化劇烈邊緣敏感，對角線方向敏感更強，(23)和(32)也對原圖像反對角線和對角線方向變化劇烈邊緣敏感，但反對角線方向敏感更強。通過增加時域窗口的大小為 7×7，9×9，11×11，13×13(由于分解圖像眾多，在此只表示5×5的情況，如圖3所示)，并分析分解重組圖像，得出結論(按照圖2分解圖像排列方式，且不包括11頻段圖像):(1)第一行表示原圖像的水平向高頻圖像，第一列表示原圖像的垂直向高頻圖像，剩余分解高頻圖像表示方向如圖4所示;(2)高頻分解方向數為2×N1-2，其中N1為時域窗口大小;(3)隨著時域窗口的增加，多元成分的干擾項越嚴重。

圖4 圖像2D-PWVD分解圖像方向表示圖

3 基于2D偽魏格納分布的SAR圖像噪聲抑制

Fukuda等［9］根據小波分解高頻系數具有不同的方向性特點，構造了3個如圖5所示的濾波窗口，來判斷高頻細節中圖像邊緣和噪聲，達到抑制噪聲的目的。

圖5 高頻子帶定向濾波器

在本文第2節闡述二維偽魏格納威爾分布分解后具有2*N-1個方向(N表示時域窗口的大小)，在Fukuda提出算法的啟示下，考慮高頻圖譜的方向性，使SAR圖像噪聲抑制后圖像邊緣特性得到保持。在本算法中，同樣采用圖5的濾波窗口，在分解的圖像的垂直向高頻圖像中，利用圖5a)濾波窗口，水平高頻圖像利用圖5b)濾波窗口，其他高頻段圖像利用圖5c)濾波窗口。具體算法如下:

(1)選擇時域窗口大小為N×N。

(2)將原圖像進行二維魏格納威爾變換，形成N×N幅分解圖像

式中:X(i，j)、Y(i，j)分別為計算后各分解圖像的各點像素值和計算前各分解圖像各像素值。

(3)針對各高頻分解圖像，計算各分解圖像的均值μ，根據如下原則

計算閾值T，將各分解圖像分成噪聲區域和可能細節區域。區分原則為:X(i，j)＜T為噪聲區域，X(i，j)≥T認為為可能細節區域。式中，X(i，j)為各分解高頻圖像的各點像素值。

(4)在噪聲區域

Xchu(i，j)=k1× X(i，j)

式中:Xchu(i，j)為處理后各點像素值;k1為噪聲抑制系數。

(5)在可能細節區域，當除中心點之外坐標中有一個點的值滿足X(i，j)≥T(圖5中深色的點)，認為當前像素值為目標邊緣細節，則 Xchu(i，j)=X(i，j)，否則認為是噪聲，Xchu(i，j)=k1×X(i，j)。

(6)將處理后的各分解圖像相加，形成SAR噪聲抑制后圖像。

4 SAR圖像相干斑抑制性能的評估

SAR圖像經過相干斑抑制算法處理后圖像質量是SAR圖像處理關心的重點，因此必須用一定的標準來衡量SAR圖像的質量。SAR圖像的應用十分廣泛，不同的用途對SAR圖像質量的要求也不完全一樣。一般來說，都要求圖像具有較好的區分臨近目標的能力和檢測目標的能力，要求圖像具有較豐富的細節信息［10］。

目前，已提出的SAR圖像斑點噪聲抑制質量定量評估指標主要包括等效視數和邊緣保持指數，下面分別介紹各個指標的定義。

1)等效視數［11］

等效視數是衡量一幅圖像相干斑噪聲相對強度的一個有效指標，等效視數越大，表明圖像上的相干斑越弱，可解譯性越好。定義為

2)邊緣保持指數［12］

邊緣保持指數定義為

式中:N為圖像的像素點數;ps，psn分別為濾波后圖像的像素值和該像素點水平或垂直方向的領域值;po，pon分別為原始圖像的像素值和該像素點水平或垂直方向的鄰域值。因此，邊緣保持指數可分為水平邊緣保持指數和垂直方向邊緣保持指數。

5 噪聲抑制試驗與結果分析

為驗證本文提出SAR圖像噪聲抑制算法的有效性，利用圣地亞哥國家實驗室提供的MiniSAR真實合成孔徑雷達數據進行試驗。

MiniSAR數據(目前公布19幅圖像)，在2005年由圣地亞哥國家實驗室獲取，采用Ka和Ku波段聚束式成像方式，最高分辨率為0.1 m×0.1 m，圖像未進行多視處理，圖像為復雜場景。

針對均值濾波算法［11］、LEE 濾波算法［13］、軟閾值小波濾波算法［1］和本文提出的噪聲抑制算法進行對比試驗，性能判別指標選用等效視數和邊緣保持系數。在試驗中，均值濾波算法采用3×3的矩形窗口，LEE濾波算法窗口大小7×7(MiniSAR數據)，軟閾值小波濾波算法中的α和β選擇為0.5，本文算法的k1=0.5，h(k，l)=1，大小為3×3的矩形窗口。由于篇幅有限，列出第9幅圖像(如圖6所示)的實驗結果(見圖7和表1)。

圖6 原始圖像

圖7為圖6圖像采用不同方法噪聲抑制后的效果圖。表1是采用不同方法的噪聲抑制性能結果。

圖7 MiniSAR數據的第9幅圖像噪聲抑制后的效果圖

表1 MiniSAR數據噪聲抑制方法性能參數表

從圖7可以看出，采用本文算法和軟閾值小波濾波算法抑制圖像從視覺效果上基本一致。

從表1可以看出，采用本文算法得到的等效視數相對于原圖像和軟閾值小波濾波算法有明顯提高，但略低于采用均值濾波和LEE濾波得到的對應參數;本文算法計算的邊緣保持系數明顯高于均值濾波和LEE濾波算法，與小波濾波算法相當。

為驗證本文算法的魯棒性，針對等效視數和邊緣保持系數指標，表2統計了本文算法相對于均值濾波算法、LEE濾波算法和軟閾值小波濾波算法對19幅MiniSAR數據實驗結果。

從表2可以看出，針對MiniSAR數據，本文算法在等效視數指標上不如均值濾波算法和LEE濾波算法優越，在邊緣保持系數指標上明顯高于其他3種算法，相對于軟閾值小波濾波算法，無論從等效視數指標還是邊緣保持系數指標，都具有一定的優越性。

從各種算法的原理上分析，本文算法將原圖像按照不同方向進行分解，并對目標細節和噪聲中可能目標細節進行分別處理，在保留圖像中盡可能多細節的同時，對噪聲進行了有效抑制。可以認為，本文算法在保持較好的邊緣保持系數的條件下，提高了SAR圖像的等效視數，達到了抑制SAR圖像噪聲的目的。

從算法效率上來分析，由于本文算法需要將原圖像分解為N幅(與時域窗口得大小相關)與原圖像同尺寸的子圖像，并需要考慮目標區域、強噪聲中可能的目標細節，因此算法效率上比均值濾波算法、LEE濾波算法和軟閾值小波濾波算法要低。

6 結束語

本文根據2D-PWVD變換的頻譜特性，提出了一種新的SAR圖像噪聲抑制算法，利用MiniSAR圖像數據進行驗證試驗，并與均值濾波、LEE濾波、軟閾值小波濾波進行對比分析，說明本文算法在抑制SAR圖像噪聲的同時，能更多保留原圖像的細節信息，證明本文算法的有效性和正確性。當然，本算法的執行效率還有待提高，以利于實際工程的具體應用。

［1］ Donoho D L，Johnstone I M.Nonlinear wavelet methods for recorvery of signal densities and spectra from indirect and noisy dat［C］//Proceedings of Symosia in Applied Mathematics.Providence Rhode Island:American Mathematical Socity，1993:173-205.

［2］ Devore R A，Lucier B J.Fast wavelet techniques for nearoptimal image processing［C］//Military Communications Conference on Fusing Command，Control and Intelligence.San Diego，CA:IEEE Press，1992:1129-1135.

［3］ Sanker B，Guler E G，Kahya Y P.Multiresolution biological transient extraction applied to respiratory crackles［J］.Computers in Biology and Medicine，1996，26(1):25-39.

［4］ 朱家兵，陶 亮，江有名，等.一種新的高分辨率SAR圖像相干斑噪聲抑制算法［J］.現代雷達，2005，27(11):54-57，74.Zhu Jiabing，Tao Liang，Jiang Youming，et al.A new speckle reduction algorithm for high resolution SAR images［J］.Modern Radar，2005 ，27(11):54-57，74.

［5］ Romberg J K，Baraniuk R G.Bayesian tree-structured image mode-ling using wavelet-domain hidden markov model［J］.IEEE Transactions on Image Processing，2001，10(7):1056-1068.

［6］ Jacobson L，Wechsler H.A paradigm for invariant object recognition of brightness，optical flow and binocular disparity images［J］.Pattern Recognition Letters，1982(1):61-68.

［7］ Jacobson L，Wechsler H.A theory for invariant object recognition in the frontoparallel plane［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1984，6(3):325-331.

［8］ Jacobson L，Wechsler H.Derivation of optical flow using a spatiotemporal-frequency approach［J］.Computer Vision，Graphics，and Image Processing，1987，38(1):29-65.

［9］ Fukuda S，Hirosawa H.Supression of speckle in synthetic aperture radar images using wavelets［J］.International Journal of Remote Sensing，1998，19(3):507-510.

［10］ 匡綱要，高 貴，蔣詠梅，等.合成孔徑雷達目標檢測理論、算法及應用［M］.長沙:國防科技大學出版社，2007.Kuang Gangyao，Gao Gui，Jiang Yongmei，et al.Synthetic aperture radar target detection theory algorithm applications［J］.Changsha:National University of Defense Technology press，2007.

［11］ Chris Olive Shaun Quegan.Understanding synthetic apearture radar images［M］.Raleigh，NC:SciTech Publishing，Inc，2004.

［12］ Sheng Yongwei，Xia Zongguo.A comprehensive evaluation of filters for radar speckle suppression［C］//International Geoscience and Remote Sensing Symposium.Lincoln，NE:IEEE Press，1996:1559-1561.

［13］ Lee J S.Digital image enhancement and noise filtering by using local statistics［J］.IEEE Transactions on Pattern and Machine Intelligence，1980，2(2):165-168.