基于矩陣對數累積量的極化合成孔徑雷達數據去噪方法

劉坤，馬文萍，劉紅英，王爽

(西安電子科技大學 智能感知與圖像理解教育部重點實驗室 智能感知與計算國際聯合研究中心，陜西西安710071)

0 引言

近些年，隨著成像雷達技術的不斷發展，越來越多的對地觀測合成孔徑雷達(SAR)系統具備了多極化或全極化成像能力。2014年4月和5月歐州太空局(ESA)與日本宇航局(JAXA)接連發射了兩顆具備全極化成像能力的新一代雷達衛星Sentinel-1A與ALOS-2/PALSAR-2，以全面提升自身的極化SAR(PolSAR)對地遙感能力。美國、加拿大等遙感強國也在按計劃逐步推進各自的PolSAR 技術發展與空、天基PolSAR 平臺的更新換代，PolSAR 領域的發展日新月異。與單極化SAR 類似，斑點噪聲抑制同樣是PolSAR 數據處理的關鍵。相干斑噪聲不僅影響PolSAR 圖像的視覺效果，更給大數據背景下Pol-SAR 圖像的智能理解、解譯帶來了困難。

經過近20 多年的不斷發展，PolSAR 去噪已經有了許多成熟的方法。如由Lee 濾波發展出來的精致極化Lee 濾波［1］、Sigma 濾波［2］、IDAN 濾波［3］等。其中精致極化Lee 濾波因其可靠穩定且對后續的分割、分類影響正面應用最為廣泛。近期隨著一些新的圖像分析方法和數學工具的發展，又涌現出許多PolSAR 去噪新方法。例如基于雙邊濾波的方法［4］等。其中在基于圖像的去噪算法中非局部方法［5］異軍突起并已逐步應用于單、多極化SAR 去噪領域。其中PPB 算法［6］利用貝葉斯定理將非局部方法應用于單極化SAR 當中，去噪效果得到本領域研究者的公認，并且在簡單假設相干矩陣僅有主對角線元素的前提下將該方法推廣到了PolSAR 去噪當中［7］。而以復Whisart 分布模型的統計分析為基礎的Pretest 算法［8］在濾波中更是利用到了相干矩陣的全部元素，作為極化濾波方法更為合理，效果也更加優秀。非局部方法作為一種空域去噪方法，同質區域選取以及相似性計算是關鍵。其相對于傳統“局部”方法能取得更好去噪效果的原因在于該方法是面向圖像塊的去噪方法，這種基于圖像塊的計算可以更好地利用圖像的結構信息，從而在圖像結構和細節信息的保持方面有突出的表現。對于PolSAR 數據來說，地物特性通常也體現為圖像上的結構和細節信息，因此該方法也適合于處理PolSAR 數據。

Anfinsen 等在復數域的梅林變換的基礎上對PolSAR 數據進行了統計分析，并在復矩陣對數累積量的基礎上提出了PolSAR 數據Goodness-of-Fit(GoF)檢測方法［9］。GoF 檢測可以通過假設檢驗來分析PolSAR 數據的兩個樣本集是否符合相同參數的Wishart 分布。但是該方法要求用來測試的樣本集包含足夠大的樣本數量，并且理論上應該是同一散射點的多次獨立觀測。Anfinsen 等［9］在模擬數據上驗證了該理論，并在真實數據上采用手工選取同質圖像塊進行了簡單的驗證，上述工作表明該統計方法適用于分析PolSAR 數據，但是還未將其真正應用于PolSAR 數據去噪或地物分類中。

對于PolSAR 數據，多視處理是常用的去噪策略。多視處理的一個重要假設就是同一數據點的各子視圖像可視為同一散射點的多次不同觀測。受此啟發本文結合非局部方法和Anfinsen 的理論［9］提出了一種基于復矩陣對數累積量的PolSAR 數據非局部去噪方法。該方法利用非局部策略選擇同質區域，然后根據同質區域數據點的子視數據集，基于復矩陣的對數累積量以GoF 方法計算同質區域的相似度，并以此為依據計算濾波權值。在非局部思想下進行的同質區域選取保留了非局部方法對圖像結構保持較好的優點，而在子視數據集基礎上計算相似度則結合了PolSAR 數據的散射特性。該方法將PolSAR 數據散射特性和圖像結構特性相結合，將基于復矩陣對數累積量的GoF 檢測方法成功應用到了PolSAR 數據處理中。該方法的優點在于同質區域選取和濾波系數的計算均是針對相干矩陣進行的，相對于僅使用主對角線上的元素或SPAN 圖像，其對極化數據的處理更具合理性。

1 基于PolSAR 復相干矩陣對數累積量的相似度計算

1.1 基于PolSAR 數據的矩陣對數累積量的計算

雷達目標的變極化效應可以用一個復二維矩陣的形式來表示，稱為極化散射矩陣或Sinclair 矩陣，用矩陣S 來表示［10］:

式中:h 和v 分別表示水平和垂直極化。散射矩陣也可用其矢量形式表示為散射矢量:

滿足完全發展的相干斑的條件時，散射矢量服從復高斯分布，其多視圖像的協方差矩陣C 為

式中:L 為視數;上標H 表示復共軛轉秩。

設d×d 復矩陣C 的概率密度函數為fC(C).

函數g(C)的梅林變換定義為

式中:Ω+表示埃爾米特正定矩陣錐。fC(C)的梅林變換定義為

當φC(s)存在時:

第v 階對數矩陣(MLM)為

第v 階矩陣對數累積量(MLC)為

式中:

MLC 可由MLM 按照下式計算得出:

PolSAR 數據中，一個數據點的n 個子視數據可視為n 個獨立同分布的復矩陣Ci，該數據點的MLM以及MLC 可以如下計算得出。

第v 階MLM:

第1 ～vp階MLM 構成vp階MLM 向量:

MLC:

將〈μv{C}〉簡記為〈μv〉，按照(9)式可得第v 階MLC 為

第1 ～vp階MLC 構成vp階MLC 向量:

1.2 基于GoF 檢測的統計檢驗值計算

Anfinsen 等基于GoF 檢測經推導得到了檢測兩個樣本集是否符合同一分布的統計檢驗值Qp的計算方式［9］:

式中:〈κ〉1與〈κ〉2分別為來自兩組樣本的MLC;K為向量

的協方差矩陣，

這里給出vp=8 時K4的表達式:記k=〈κ〉2，則

式中:

2 結合非局部均值的去噪算法

2.1 基于非局部均值的同質區域選取

非局部均值方法通過兩個圖像塊之間所有對應數據點的相似程度總和來衡量圖像塊中心像素的相似程度。圖像中每一個像素點的濾波結果由搜索窗內同質區域中所有像素點按照與待濾波點的相似程度加權得到。本文采用Pretest 方法［3］中的閾值方法來選取同質區域。為了保證濾波結果能夠保持圖像的結構特征，同質區域選取是在空域多視平均過后的圖像上進行的。任意兩個圖像塊之間可計算統計檢測值:

式中:d 為相干矩陣的階數;w 為一個圖像塊中的點數;Xi與Yi分別為兩個圖像塊對應位置處的相干矩陣。如圖1所示。

圖1 非局部方法選取同質區域示意圖Fig.1 Homogeneous region selected using non-local method

閾值可以通過下式進行估計:

對搜索窗內所有非待濾波點計算ln H，當

時將該點記入待濾波點的同質區域參與濾波。

2.2 計算同質點與待濾波點的相似度并加權去噪

非局部方法中兩點之間的相似度也是通過兩點所在圖像塊之間的相似度進行計算的。而對于未經過多視處理的PolSAR 數據而言，兩個點對應的兩組子視數據集同樣可以很好地反應出相似程度。

包含n 個子視數據的數據點的MLM 可利用其各子視C 矩陣Ci(i=1，…，n)按照(10)式的計算方法計算得出。這樣其MLC 可按照(12)式的計算方法計算得出。數據點1 和2 的相似程度可以由Qp反映出來，Qp按照(14)式計算得出。

濾波權值為

式中:α ＞0 為常數，用于調整濾波效果。

本文算法步驟如下:

步驟1 對全圖數據求所需要的行列式等后續計算需要的變量。

1)計算空域多視平均并選擇K 矩陣的階數vp.本文中選用vp=4.

2)計算子視數據以及多視數據的C 矩陣的行列式。

3)由子視數據行列式按照(10)式計算2vp階MLM 向量，按照(12)式和(13)式計算2vp階MLC向量并按照(16)式計算vp階K 矩陣。本文中使用(17)式計算4 階K 矩陣。

步驟2 對一個待濾波點進行濾波。

1)對待濾波點選擇同質區域。按照(23)式選取合理的閾值ln Ht. 由空域多視平均結果在搜索窗內按照(22)式計算搜索窗內所有點與待濾波點的統計檢驗值ln H，并按照(24)式選擇該待濾波點的同質區域。

2)計算同質點的濾波權值。由同質區域中所有點的MLC 向量以及待濾波點處的K 矩陣按照(14)式計算Qp，并按照(25)式計算濾波權值ω.

3)對待濾波點加權濾波。該待濾波點C 矩陣濾波結果為

式中:m 為同質區域中的數據點數。

步驟3 循環下一待濾波點，按照步驟2 進行濾波。

本算法的運算復雜性主要由Pretest 方法選取同質區域的運算復雜性決定，在文獻［3］中有該方法的快速算法，按照該快速算法如將(22)式中的計算視為一步，M×N 的圖像在窗口大小為f 時Pretest去噪方法的計算規模可降低至O(M×N×f2). 采用該快速方法后本文算法的時間復雜度與此相同。

3 去噪效果對比分析

大多數成像雷達的方位分辨率高于距離分辨率，因此盡管兩個方向都可進行多視處理，但一般都選在方位向。目前學術研究中常用的公開PolSAR數據有很多，大多數都是在頻域經過多視處理之后的多視復數據。加拿大的機載全極化SAR 系統Convair-580 公開發布的數據距離分辨率約為4 m，方位分辨率約為0.4 m，在對其進行處理前通常先在方位向進行空域的10 視平均，較合適用來測試本算法的實際效果并與其他方法進行對比。

這里使用兩組真實的Convair-580 獲取的數據進行實驗，其中第一組Ottawa 數據是2001年6月26日獲得的渥太華地區數據(見圖2)，多視處理前數據大小為222×3 429(距離向×方位向)，取其中222×3 420 大小的數據作為實驗數據，10 視處理后大小為222×342. 另一組NRCAN 數據來源于加拿大自然資源部(NRCAN)網站(見圖3)，因其數據頭不全，僅知數據獲取于2000年6月20日，具體地點位置并不清楚。其數據大小為1 024×11 000(距離向×方位向)，本文截取其中256×2 560(距離向×方位向)大小的數據作為實驗之用，10 視處理后大小為256×256. 選取的區域中有市區、道路、不同類型植被，并且有成對的電線塔架作為大目標以及道路上的車輛作為點目標，十分適于用來主觀評價去噪效果好壞。

圖2 Ottawa 數據去噪效果對比Fig.2 Comparison of different filters on Ottawa dataset

圖3 NRCAN 數據去噪效果對比Fig.3 Comparison of different filters on NRCAN dataset

等效視數(ENL)是評價SAR 數據區域平滑程度的常用指標，其定義為標準差與均值之比。ENL數值越大表示用于計算ENL 區域的平滑效果越好。本文在去噪后的SPAN 數據上計算ENL，并在計算ENL 時為了實驗結果更加客觀，選擇了不同位置的多個同質區塊。

在評價去噪算法時只考慮同質區域的平滑效果是不夠的，通常在追求同質區域更加平滑、受相干斑影響更小的同時還需要圖像的結構信息、細節等有更小的失真。2012年Mittal 等提出了無參考圖像的空域評價指標(BRISQUE)［11］，文獻［12］將該指標用于PolSAR 圖像相干斑抑制效果評價。BRISQUE 可以在空域評價圖像的質量，該指標是通過計算圖像局部正規化亮度系數的統計信息來評價圖像的失真度的，其優點在于無需參考圖像，非常適合像SAR 斑點噪聲抑制這樣沒有真實參考圖像的問題。因此，本文選擇該指標對Pauli 偽彩圖進行評價。BRISQUE 的數值在0 ～100 之間，結果越小表示圖像的質量越好。BRISQUE 的計算公式如下:

式中:ξ 表示尺度;σ 表示方差;Γ(·)表示Gamma 函數;ψ=σ

本文實驗在Intel Core i3 CPU(3.2 GHz)，4G RAM 的環境下進行。實驗所用軟件為Matlab R2013a 以及PolSARPro 4.2.0. 對比算法包括均值方法(Boxcar)，窗口大小為5;精致極化Lee 濾波［1］，窗口大小為7;Sigma 濾波［2］，Sigma 值為0.9，目標窗口為3，濾波窗口為9;Pretest［8］濾波，窗口大小為3，搜索窗大小為9，閾值ln Ht為-120. 前面3 組對比算法均使用歐洲太空局的PolSAR Pro 軟件實現。本文算法中非局部窗口大小為3，搜索窗大小為9，閾值ln Ht為-120，濾波權值計算式中系數α 為2.

兩組數據的去噪結果如圖2和圖3所示。在每組數據的原圖中標出了用于測試ENL 值的圖塊A和圖塊B. 各種去噪算法的噪聲抑制性能以ENL 和BRISQUE 作為標準整理在表1中。

表1 噪聲抑制效果對比Tab.1 Comparison of speckle filtering effects of different methods in terms of ENL and BRISQUE

從對比實驗可以看出，本文方法在抑制噪聲的同時，很好地保持了圖像的結構和細節信息。Sigma濾波、Pretest 濾波以及本文方法對同質區域的平滑效果較好，總體大致相當。本文方法與Pretest 方法相比具有相當的同質區域平滑能力，并且圖像失真更小。在圖像比較復雜的區域，如NRCAN 右下部分以及Ottawa 數據左上部分的城區(住宅區)，都可直觀看出本文方法在細節和結構信息保持方面較Sigma 濾波、Pretest 濾波更具優勢，其中Sigma 濾波存在“過平滑”的問題。這一點在BRISQUE 值上得到了體現。在兩個數據集上本文方法去噪結果的BRISQUE 值比其他方法都要低出許多。而ENL 所體現出來的同質區域的平滑效果則說明了本文方法中相似度計算以及濾波權值計算的合理性。

本文方法對目標特性的保持也較好，例如NRCAN 圖像中由左中到右上排成一列的電線塔架區域，由濾波結果可以看出本文方法的去噪結果相比原始圖像的失真最小。而觀察圖中的小目標，例如NRCAN 數據中央部分相互交叉的兩條道路上的車輛等，本文去噪結果中這些點目標也得到了很好的保持。另外本文方法對兩個數據集中道路邊緣的信息保持得也很好，邊緣明確并且道路并未被膨脹或腐蝕。以上這些都表明本文方法能夠在有效抑制噪聲的同時，可以更好地保持圖像的結構信息以及細節信息。

PolSAR 去噪算法的另一個關鍵指標就是對極化相關性保持的能力。為了檢測本文算法是否能保持數據的極化相關性，在NRCAN 數據中選取了A、B、C 3 個區域，分別位于植被、電線塔架以及居民區，代表3類典型地物(如圖4，因區域較小故用紅色實心方塊標示)。在原始數據以及本文方法去噪結果數據中分別計算這3 個區域散射矩陣的平均值，并繪制同極化、交叉極化特征圖，如圖5和圖6所示。圖5、圖6中相同子序號的圖像對應同一圖塊去噪前后的對應同極化或交叉極化特征。對比可見，本文算法去噪前后對各類地物都能很好地保持極化相關性。

圖4 用于測試極化特征的圖塊位置Fig.4 Patches used to test the polarization signature

圖5 去噪前的同極化與交叉極化特征圖Fig.5 The polarization signatures before filtering by the proposed method

4 結論

PolSAR 數據相干矩陣對數累積量的計算方法自2011年提出以來一直是PolSAR 數據處理領域的一個新亮點，但理論研究較多，而在噪聲抑制、分割、分類、目標檢測識別等應用方向卻少有突破。本文將該理論和非局部去噪思想相結合，提出了一種新的PolSAR 去噪算法。新方法使用了相干矩陣的行列式來選擇同質區域以及計算相似度，相對于僅使用主對角線上的元素或SPAN 圖像更具合理性。本文提出的方法首次直接將復矩陣對數累積量應用于PolSAR 斑點噪聲抑制中，文中從噪聲平滑、圖像保真、極化信息保持三方面對比了新方法與其他Pol-SAR 去噪方法的性能，展現了該方法的特性與優點。本文方法為將PolSAR 數據相干矩陣對數累積量的計算理論進一步應用到PolSAR 數據處理的更多應用提供了參考。

圖6 去噪后的同極化與交叉極化特征圖Fig.6 The polarization signatures after filtering by the proposed method

References)

［1］ Lee J S，Grunes M R，Grandi G. Polarimetric SARspeckle filtering and its implication for classification［J］.IEEE Transactions on Geoscience and Remote Sensing，1999，37(5):2363 -2373.

［2］ Lee J S，Wen J H，Ainsworth T L，et al. Improvedsigma filter for speckle filtering of SAR imagery［J］.IEEE Transactions on Geoscience and Remote Sensing，2009，47(1):202 -213.

［3］ Vasile G，Trouvé E，Lee J S，et al. Intensity-driven adaptiveneighborhood technique for polarimetric and interferometric SAR parameters estimation［J］. IEEE Transactions on Geoscience and Remote Sensing，2006，44(6):1609 -1621.

［4］ 王爽，于佳平，劉坤. 基于雙邊濾波的極化SAR 相干斑抑制［J］.雷達學報，2014，3(1):35 -44.WANG Shuang，YU Jia-ping，LIU Kun. Polarimetric SAR speckle reduction based on bilateral filtering［J］.Journal of Radars，2014，3(1):35 -44.(in Chinese)

［5］ Buades A，Coll B，Morel J M. A non-local algorithm for image denoising［C］∥CVPR 2005. San Diego，CA:IEEE，2005:60-65.

［6］ Deledalle C A，Denis L，Tupin F. Iterative weighted maximum likelihood denoising with probabilistic patch-based weights［J］.IEEE Transactions on Geoscience and Remote Sensing，2009，18(12):2661 -2672.

［7］ Deledalle C A，Tupin F，Denis L. Polarimetric SAR estimation based on non-local means［C］∥IGARSS 2010. Honolulu，HI:IEEE，2010:2515 -2518.

［8］ Jiong C，Yilun C，Wentao A. Nonlocal filtering for polarimetric SAR data:a pretest approach［J］. IEEE Transactions on Geoscience and Remote Sensing，2011，49(5):1744 -1754.

［9］ Anfinsen S N，Doulgeris A P，Eltoft T. Goodness-of-Fit tests for multilook polarimetric radar data based on the mellin transform［J］. IEEE Transactions on Geoscience and Remote Sensing，2011，49(7):2764 -2781.

［10］ Lee J S. Polarimetric radar imaging［M］. New York:Taylor ＆Francis Group，2009:63 -68.

［11］ Mittal A，Moorthy A K，Bovik A C. No-reference image quality assessment in the spatial domain［J］.IEEE Transactions on Image Processing，2012，21(12):4695 -4708.

［12］ Torres L，Sant’Anna S J S，da Costa Freitas C，et al. Specklereduction in polarimetric SAR imagery with stochastic distances and nonlocal means［J］.Pattern Recognition，2014，47(1):141-157.