劉 聰，李言俊，劉保元

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

0 引言

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

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

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

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

1 二維魏格納分布理論

Wigner于1932年首先提出了Wigner分布的概念，并應(yīng)用于量子力學(xué)領(lǐng)域。1948年Ville將Wigner分布引入信號(hào)分析領(lǐng)域，應(yīng)用于信號(hào)的檢測和信號(hào)細(xì)節(jié)分類。因此Wigner分布又稱為Wigner-Ville分布(WVD)，Wigner-Ville分布屬于雙線性時(shí)頻分析Cohen＇s類的特殊形式，具體表達(dá)式為

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

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

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

離散形式定義為

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

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

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

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

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

圖1 圖像數(shù)據(jù)的2D-PWVD分解流程

圖2 分解圖像的排列方式

選取時(shí)域窗口函數(shù)h(k，l)=1(矩形窗)，大小為5×5，通過圖1的變換重組流程，將lenna原始圖像分解為25個(gè)不同的頻段圖(圖3)。從分解重組圖可以看出，圖3(11)相對于lenna原始圖像，圖像邊緣模糊，可以理解為是原始圖像經(jīng)過2D-PWVD變換重組的低通濾波圖像，分解圖像(12)～(55)是原始圖像邊緣的信息體現(xiàn)，可以理解為是原始圖像經(jīng)過2D-PWVD變換重組的高通濾波圖像。仔細(xì)分析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圖像經(jīng)過2D-PWVD變換重組，實(shí)際上形成了13個(gè)分解圖像。推廣到所有時(shí)域窗口的情況下，分解重組后形成的圖像數(shù)為(N1×N2+1)/2。

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

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

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

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

Fukuda等［9］根據(jù)小波分解高頻系數(shù)具有不同的方向性特點(diǎn)，構(gòu)造了3個(gè)如圖5所示的濾波窗口，來判斷高頻細(xì)節(jié)中圖像邊緣和噪聲，達(dá)到抑制噪聲的目的。

圖5 高頻子帶定向?yàn)V波器

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

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

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

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

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

計(jì)算閾值T，將各分解圖像分成噪聲區(qū)域和可能細(xì)節(jié)區(qū)域。區(qū)分原則為:X(i，j)＜T為噪聲區(qū)域，X(i，j)≥T認(rèn)為為可能細(xì)節(jié)區(qū)域。式中，X(i，j)為各分解高頻圖像的各點(diǎn)像素值。

(4)在噪聲區(qū)域

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

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

(5)在可能細(xì)節(jié)區(qū)域，當(dāng)除中心點(diǎn)之外坐標(biāo)中有一個(gè)點(diǎn)的值滿足X(i，j)≥T(圖5中深色的點(diǎn))，認(rèn)為當(dāng)前像素值為目標(biāo)邊緣細(xì)節(jié)，則 Xchu(i，j)=X(i，j)，否則認(rèn)為是噪聲，Xchu(i，j)=k1×X(i，j)。

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

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

SAR圖像經(jīng)過相干斑抑制算法處理后圖像質(zhì)量是SAR圖像處理關(guān)心的重點(diǎn)，因此必須用一定的標(biāo)準(zhǔn)來衡量SAR圖像的質(zhì)量。SAR圖像的應(yīng)用十分廣泛，不同的用途對SAR圖像質(zhì)量的要求也不完全一樣。一般來說，都要求圖像具有較好的區(qū)分臨近目標(biāo)的能力和檢測目標(biāo)的能力，要求圖像具有較豐富的細(xì)節(jié)信息［10］。

目前，已提出的SAR圖像斑點(diǎn)噪聲抑制質(zhì)量定量評估指標(biāo)主要包括等效視數(shù)和邊緣保持指數(shù)，下面分別介紹各個(gè)指標(biāo)的定義。

1)等效視數(shù)［11］

等效視數(shù)是衡量一幅圖像相干斑噪聲相對強(qiáng)度的一個(gè)有效指標(biāo)，等效視數(shù)越大，表明圖像上的相干斑越弱，可解譯性越好。定義為

2)邊緣保持指數(shù)［12］

邊緣保持指數(shù)定義為

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

5 噪聲抑制試驗(yàn)與結(jié)果分析

為驗(yàn)證本文提出SAR圖像噪聲抑制算法的有效性，利用圣地亞哥國家實(shí)驗(yàn)室提供的MiniSAR真實(shí)合成孔徑雷達(dá)數(shù)據(jù)進(jìn)行試驗(yàn)。

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

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

圖6 原始圖像

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

圖7 MiniSAR數(shù)據(jù)的第9幅圖像噪聲抑制后的效果圖

表1 MiniSAR數(shù)據(jù)噪聲抑制方法性能參數(shù)表

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

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

為驗(yàn)證本文算法的魯棒性，針對等效視數(shù)和邊緣保持系數(shù)指標(biāo)，表2統(tǒng)計(jì)了本文算法相對于均值濾波算法、LEE濾波算法和軟閾值小波濾波算法對19幅MiniSAR數(shù)據(jù)實(shí)驗(yàn)結(jié)果。

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

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

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

6 結(jié)束語

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

［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］.現(xiàn)代雷達(dá)，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］ 匡綱要，高 貴，蔣詠梅，等.合成孔徑雷達(dá)目標(biāo)檢測理論、算法及應(yīng)用［M］.長沙:國防科技大學(xué)出版社，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.