改進分塊局部最佳維納濾波算法的干涉相位濾波*

汪 洋，黃海風，董 臻，吳曼青,2

改進分塊局部最佳維納濾波算法的干涉相位濾波*

汪 洋1，黃海風1，董 臻1，吳曼青1,2

(1.國防科技大學 電子科學與工程學院， 湖南 長沙 410073；

2.中國電子科技集團， 北京 100000)

針對合成孔徑雷達干涉相位濾波問題，提出了一種改進的分塊局部最佳維納濾波算法。該算法是加性高斯白噪聲下的線性最小均方誤差估計，利用目前圖像濾波最前沿的技術——非局部技術，來聯合估計圖像的一、二階矩。針對干涉相位中噪聲的空變性，在應用中提出了兩點改進：估計噪聲的標準差時，用均值代替中值；根據噪聲標準差的最大值和均值的比值，自適應地確定類的數量。仿真和實測數據表明，改進后的分塊局部最佳維納濾波算法是有效的，并優于其他三種算法。

干涉相位濾波；分塊局部最佳維納濾波；線性最小均方誤差估計

(1.CollegeofElectronicScienceandEngineering,NationalUniversityofDefenseTechnology,Changsha410073,China;

2.ChinaElectronicsTechnologyGroupCorporation,Beijing100000,China)

合成孔徑雷達(SyntheticApertureRadar,SAR)由于其全天時、全天候、主動遙感的優良特性，在軍事偵察、國民經濟建設和科學研究中具有廣泛的應用，干涉測量[1]便是其中之一。干涉相位是合成孔徑雷達干涉測量(SyntheticApertureRadarInterferometry,InSAR)中最重要的物理量，其質量的好壞將決定最終產品——數字高程模型(DigitalElevationModel,DEM)的精度。然而，受到相關因素[2-3]的影響(時間、空間、體散射、噪聲等)，干涉相位中總是存在嚴重的噪聲。這些噪聲不僅會引入殘差點，還會破壞干涉條紋的分布增加了相位解纏的難度，最終導致DEM精度的降低，因此必須予以濾除。

總的來說，干涉相位的濾波方法可以歸結為兩類：相位域和變換域。相位域的方法對干涉相位直接濾波而不作變換，代表算法如Lee濾波[4]及其改進算法[5]、旋濾波[6]等。變換域的方法是將干涉相位變換到其他域濾波，代表算法如Goldstein濾波[7]及其改進算法[8]、小波算法[9]、小波包算法[10]。這兩大類算法都試圖根據信號與噪聲不同的統計特性將其區分開，達到濾除噪聲保持信號的目的。Lee濾波[4]是一種各向異性的濾波算法，它利用局部統計特性和自適應窗口來濾波；旋濾波[6]根據條紋與噪聲在條紋法線方向和切線方向的統計特性，采用自適應窗口濾波；Goldstein[7]最先將頻域的方法引入干涉相位濾波，算法把干涉相位圖從相位域轉換到頻域，對頻譜進行平滑。由于噪聲在頻域是寬帶信號，而有用信號是窄帶信號，對圖像進行頻域的低通濾波就可以實現去噪。然而，算法中的濾波參數卻是固定的，針對這一點Baran等[8]引入干涉相干系數，使得算法能根據干涉相干系數的大小進行自適應濾波。針對低相干區域的濾波問題，Lopez-Martinez等[9]提出了小波變換的方法。該方法能夠分解圖像的頻譜，有望將噪聲和高頻信號進行一定的區分。小波包變換能夠將圖像的高頻部分進一步細分，根據這一性質Zha等[10]提出了基于小波包變換和維納濾波的干涉相位濾波方法。

2005年，Buades等[11]提出了一種非局部平均濾波(Non-LocalMeans,NLM)的算法，開創了基于“非局部”濾波的先河。非局部思想認為，圖像中存在大量相似的小塊(這些小塊在空間上可以是不相鄰的，所以稱之為非局部)，將這些小塊基于某種相似性準則進行聚類后聯合濾波，可以提升傳統局部濾波算法的性能。目前已有許多學者對該方法進行了改進與拓展，這種非局部去噪的思想在圖像(視頻)處理[12-14]、醫學影像[15-16]等多個領域都已得到應用。根據這一思想，Chatterjee等[17]提出了加性高斯白噪聲模型下的線性最小均方誤差估計(LinearMinimumMeanSquareEstimation,LMMSE)——分塊局部最佳維納濾波算法(Patch-basedLocallyOptimalWiener,PLOW)，并推導了該算法和非局部平均濾波算法的關系。本文將PLOW算法應用于干涉相位濾波，并根據干涉相位中噪聲空變的特點進行了改進，使得算法能自適應地抑制噪聲的同時保持條紋細節。

1 相位模型

1.1 相位域模型

干涉相位是對配準后的兩幅SAR圖像的共軛乘積取相位得到的，其質量受兩幅SAR圖像之間相關系數大小的影響：

(1)

許多學者基于高斯散射模型，推導出了多視情況下分布式目標干涉相位的概率密度函數[18]：

(2)

(3)

單視情況下(n=1)，式(3)積分可以簡化為[1]：

(4)

其中Li2(·)是歐拉以2為底的對數：

(5)

基于式(2)，Lee等[4]推導出了干涉相位的加性噪聲模型：

θz=θx+ν

(6)

其中，θz是干涉相位的測量值，θx是不含噪聲的干涉相位，ν是0均值噪聲并與θx獨立。

1.2 復數域模型

由于干涉相位是-π到π的周期分布，如果直接在相位域濾波會消除相位的跳變點。然而，相位跳變的地方往往是信號的高頻部分，應該對其保留以便進行正確的相位解纏。為了解決這個問題，可以將干涉相位變換到復數域處理，其復數域表達式為：

ejθz=cos(θz)+jsin(θz)

(7)

結合式(6)，Lopez-Martinez等推導了復數域干涉相位實部和虛部的表達式分別為[9]：

cos(θz)=Nccos(θx)+νr

(8)

sin(θz)=Ncsin(θx)+νi

(9)

2 分塊局部最佳維納濾波算法

2.1 算法原理

NLM[11]算法充分利用了圖像中的自相似性，即在圖像中往往會存在一些空間不相鄰但彼此非常相似的圖像塊。接著，通過一個簡單的加權平均來估計一個像素，權值唯一地依賴于兩個圖像塊的相似性而與位置無關。與傳統局部濾波方法比較，NLM能充分利用圖像的冗余性，在有效抑制噪聲的同時很好地保留圖像的紋理結構。干涉相位的特點是具有大量周期性重復出現的條紋，圖像冗余性多。從最大似然估計的角度上講，通過增加樣本數可以降低估計的方差以達到改善估計的性能。因此，將非局部的思想應用到干涉相位濾波中是合適和有效的。

在追求濾波極限性能的驅使下[19-20]，Chatterjee等[17]提出了PLOW算法。相比NLM算法，PLOW算法有兩點主要的不同。第一，NLM是以像素點為濾波單位，PLOW則是以圖像塊為單位。第二，NLM只利用了光相似信息(photometricsimilarity)，PLOW不僅利用了光相似信息還用到了幾何相似信息(geometricsimilarity)。在PLOW中，基于單個像素點的加性噪聲模型可以寫成基于圖像塊的形式：

yi=zi+ηi

(10)

其中，I(i∈I)表示圖像分塊集，zi表示無噪圖像塊，ηi是加性噪聲，yi是含噪圖像塊。利用克拉美勞界限推導出基于式(10)的濾波性能極限為：

(11)

第一，圖像中包含幾何相似的塊，假定這些塊都具有相同的分布函數，因此應該利用基于特征的方法將含有不同幾何信息的圖像塊聚類到不同的類中去。正確的聚類是算法的基礎，因此特征必須對圖像對比度、噪聲等魯棒性進行特征提取。完成特征提取后，采用K均值算法進行聚類，而類的數量也將影響算法的性能。

第二，完成聚類后，用維納濾波器進行濾波：

(12)

(13)

第三，為了避免圖像塊之間的不連續，分塊會有重疊。因此對于單個像素點會有多個估計值，文獻[17]采用一種基于線性最小均方誤差估計的方法對所有估計值進行了加權融合處理，從而得到最終的估計值。

2.2 針對干涉濾波的改進

2.2.1 自適應噪聲方差估計

干涉相位中的噪聲是空變的，而文章中的方法是基于非空變噪聲的估計，為此文獻[21]將中值估計的結果乘以一個因子以體現噪聲的空變性，即：

(14)

(15)

其中，mean(·)是取均值，Y是圖像Y的梯度圖像[17]。對比式(14)和式(15)可以看出，式(14)實際是對空變噪聲方差的“過估計”,式(15)也是基于這樣一種思想，但具有自適應性。

2.2.2 自適應類數量估計

基于估計的相關原理，當類的數量K取的過少時，幾何上不相似的塊被分到一起，導致基于類的參數估計不準確；而當K取的過多時，每一個類中的塊較少，導致魯棒參數估計不穩健。文獻[17]中取K=15，并發現濾波的結果對K取值的變化不是很敏感。

實驗中發現，K固定地取15對所有圖像都是最優的。其次，K的取值對于K均值聚類的結果有很大的影響，可能導致算法陷入局部而非全局最優。對于結構復雜的圖像，需要較多的類以便對圖像分塊進行合理正確的聚類；而對于結構簡單的圖像，只需較少的類便可完成聚類，如果仍用較多的類，算法將會對圖像進行更精細的聚類，類與類之間的距離閾值會變小。當小到一定程度的時候，聚類會由于噪聲的影響變得不穩健。

干涉相位是(-π,π]的周期分布，結構簡單；干涉相位中存在大量噪聲和地形起伏使得結構復雜化。結合干涉相位的特點、K的取值原則，我們給出K的新的取值方法：

(16)

式(16)中，用最大值除以均值實質是對噪聲空變性的定量估計：當噪聲變化大時，對分塊相似性的影響較大，K應該取較大的值；當噪聲變化較小時，對分塊相似性的影響較小，K應該取較小的值。因此，新的方法兼顧了干涉相位取值范圍、噪聲的空變性等因素，并且能夠根據不同的數據自適應地取值。

3 實驗結果及分析

本節將從仿真和實測數據驗證算法的有效性，并將結果和經典的Goldstein算法[7]、Lee算法[4]和PLOW算法[17]作對比，驗證算法的優越性。在定量評價指標方面，采用經典的均方誤差估計(MeanSquareEstimation,MSE)(越小越好)、殘差點數量(NumberOfResidues,NOR)[22](越小越好)、邊緣保持指數(EdgePreservingIndex,EPI)[23](越接近1越好)、結構相似性(StructuralSimilarity，SSIM)[24]矩陣(越大越好)均值和算法運行時間。SSIM陣利用人類視覺系統的特性，克服了傳統MSE、峰值信噪比(PeakSignal-to-NoiseRatio,PSNR)等經典評價指標與視覺感知質量不相符的問題，著重描述兩幅圖像的結構相似性大小而非逐點相似性大小。因此，SSIM可以衡量算法的細節保持能力，詳細信息請參閱文獻[24]。

3.1 仿真數據

不含噪聲的真實相位如圖1(a)所示，兩個螺旋[22]纏繞在一起，螺旋的邊緣處像素值發生了劇烈跳變，圖像大小為257×257像素。為了模擬干涉相位中的空變噪聲，簡單地將圖像平均分成四個部分，分別添加標準差為0.3，0.5，0.7和0.9弧度的高斯白噪聲，經2π纏繞后生成的干涉相位如圖1(b)所示。從中可以清楚地看到，受噪聲的影響邊緣變得模糊。

(a)無噪相位(a) Real phase

(b)含噪相位(b) Interferometric phase圖1 真實相位和干涉相位Fig.1 Real phase and interferometric phase

由Goldstein算法、Lee算法、PLOW算法和本文算法得到的濾波結果分別如圖2(a)、(b)、(c)和(d)所示。從圖中可以看到，由Goldstein算法得到的結果中還存在明顯的噪聲，從圖1(a)的下半部分可以看出。相比無噪相位，垂直方向條紋邊緣上的噪聲依然沒有得到很好的平滑，說明Goldstein算法是欠濾波的。Lee算法得到的結果要優于Goldstein算法，從圖2(b)中可以看到，圖像下半部分的噪聲相比Goldstein算法得到了較好的抑制，而這部分噪聲在PLOW算法(圖2(c))中得到了更進一步地抑制。然而，圖2(b)、(c)中垂直向條紋上的噪聲仍然較明顯。本文算法得到的結果如圖2(d)所示，不但垂直向條紋上的噪聲得到了較好的平滑，而且圖像的細節信息也得到了較好的保留，其結果是4種算法中最優的。定量的評價指標如表1所示。在MSE方面，Goldstein算法最差，這是由于欠 濾波導致的，Lee算法明顯優于Goldstein算法，而本文提出的算法是最優的。在NOR方面，無噪干涉相位有48個，含噪相位有1086個，除了Goldstein算法外，Lee、PLOW和本文算法差別不大。在邊緣保持能力上本文算法是最優的，其EPI值最接近1，PLOW算法優于Lee算法排在第二，Goldstein算法的邊緣保持能力最差。圖像的結構信息是圖像最重要的信息，如果保持得當，原始的圖像信息幾乎可以通過一個簡單的線性逆變換恢復出來。因此，任何圖像處理技術都要盡量避免破壞結構信息。那么，從平均結構相似度(MeanStructuralSimilarity,MSSIM)上說，本文算法優于其他三種算法，對于圖像結構的保護是最好的，這對相位解纏和DEM反演是極其有利的。最后，從算法的用時上看，本文算法和最快的Goldstein算法相近，都明顯優于Lee算法和PLOW算法。

(a) Goldstein算法得到的濾波結果(a) Filtered phase by Goldstein filter

(b) Lee算法得到的濾波結果(b) Filtered phase by Lee filter

(c) PLOW算法得到的濾波結果(c) Filtered phase by PLOW filter

(d) 本文算法得到的濾波結果(d) Filtered phase by the proposed filter圖2 不同算法的識別結果Fig.2 Filtered phases by different methods

算法MSENOREPIMSSIM用時/s無濾波0．641810865．91680．2704—Goldstein0．46351062．25610．68039．2Lee0．2140601．43720．765659．9PLOW0．1739501．31150．813648．7本文算法0．1546461．14580．856410．1

3.2 實測數據

本數據的成像區域是意大利的Etna火山，截取大小為100×100像素的干涉相位，如圖3(a)所示。干涉條紋很密集并且噪聲影響嚴重，這對考驗算法在低相干、密集條紋區域的性能具有重要意義。

由Goldstein算法、Lee算法、PLOW算法和本文算法得到的濾波結果分別如圖3(b)、(c)、(d)和(e)所示。對比結果可以看出：圖3(b)中的噪聲有一定程度的濾除，條紋的清晰度有一定程度的改善，但噪聲的影響依然嚴重。而從Lee算法得到的結果可以看出，由于密集條紋的影響，條紋發生了融合和斷裂，并且無法得到有效的干涉條紋。圖3(d)中的噪聲較圖3(b)得到了進一步地抑制，條紋的清晰度較圖3(b)有了改善。而從本文算法得到的濾波結果(圖3(e))中可以看出，條紋清晰可辨并且未發生斷裂或者融合，說明本文算法在抑制噪聲的同時能夠較好地保持圖像的細節。

(a)干涉相位(a) Interferometric phase

(b) Goldstein算法得到的濾波結果(b) Filtered phase by Goldstein filter

(c)Lee算法得到的濾波結果(c)Filtered phase by Lee filter

(d)PLOW算法得到的濾波結果(d) Filtered phase by PLOW filter

(e) 本文算法得到的濾波結果(e) Filtered phase by the proposed filter圖3 干涉相位和不同算法得到的濾波結果Fig.3 Interferometric phases and filteredones obtained by different algorithms

算法NORNOR減少百分比/%用時/s無濾波20610—Goldstein129237．30．9Lee75663．311．2PLOW71365．48．6本文算法12893．81．2

在定量評價方面，由于沒有真值相位，我們只能用NOR和算法運行時間來評估，其結果如表2所示。原始干涉相位含有2061個殘差點，Goldstein算法在不破壞條紋的前提下只濾除了小部分噪聲，而Lee算法雖然將殘差點的數量降低了63.3%，但這是以破壞干涉條紋為代價的，說明Lee算法的自適應窗口并不能對密集條紋進行有效的擬合，因而不適合密集條紋的濾波。由于采用了非局部的思想，PLOW算法無論從視覺效果上還是殘差點抑制上都要優于Goldstein和Lee算法，而本文提出的算法將殘差點數量降低了93.8%，并且在運行時間上較PLOW算法有了很大提高，接近最快的Goldstein算法，這說明本文提出的算法是快速有效的。

4 結論

基于分塊聯合濾波的非局部平均思想是目前濾波技術的研究熱點，這一類算法利用圖像中的冗余信息進行聯合濾波，本文將這一方法應用到合成孔徑雷達干涉相位濾波中，提出了改進的PLOW算法。改進后的算法能在空變噪聲的影響下，自適應地估計噪聲的方差和類的數量。相比原算法，改進后的算法在速度和精度上都有提高。仿真和實測數據的結果表明，本文算法不僅在濾除噪聲的同時保持圖像細節的綜合能力上優于Goldstein算法、Lee算法和PLOW算法，而且能較好地保持圖像的結構信息，這對后續的干涉處理是很有利的。

由于沒有真值，對實測數據處理結果的定量評估只能用NOR，從本文的分析看出，NOR只能作為一種粗略的評價準則。基于此，本文的下一步工作是將算法應用于更多的實測數據處理，同時研究更加精準的濾波評價準則。

References)

[1]BamlerR,HartlP.Syntheticapertureradarinterferometry[J].InverseProblem, 1998, 14(4):R1-R54.

[2] 王超, 張紅, 劉智, 等. 星載合成孔徑雷達干涉[M]. 北京：科學出版社，2002.

WANGChao,ZHANGHong,LIUZhi,etal.Spacebornesyntheticapertureradarinterferometry[M].Beijing:SciencePublishingHouse, 2002. (inChinese)

[3]ZebkerHA,VillasenorJ.Decorrelationininterferometricradarechoes[J].IEEETransactionsonGeoscienceandRemoteSensing, 1992, 30(5):950-959.

[4]LeeJS,PapathanassiouKP.AnewtechniquefornoisefilteringofSARinterferometricphaseimages[J].IEEETransactionsonGeoscienceandRemoteSensing, 1998, 36(5):1456-1465.

[5]ChaoCF,ChenKS,LeeJS.RefinedfilteringofinterferometricphasefromInSARdata[J].IEEETransactionsonGeoscienceandRemoteSensing, 2013, 51(12):5315-5323.

[6]YuQF,YangX,FuSH,etal.Anadaptivecontouredwindowfilterforinterferometricsyntheticapertureradar[J].IEEEGeoscienceandRemoteSensingLetters, 2007, 4(1):23-26.

[7]GoldsteinRM,WernerCL.Radarinterferogramfilteringforgeophysicalapplications[J].GeophysicalResearchLetters, 1998, 25(21):4035-4038.

[8]BaranI,StewartMP,KampesB,etal.AmodificationtotheGoldsteinradarinterferogramfilter[J].IEEETransactionsonGeoscienceandRemoteSensing, 2003, 41(9):2114-2118.

[9]MartinezCL,FabregasX.ModelingandreductionofSARinterferometricphasenoiseinthewaveletdomain[J].IEEETransactionsonGeoscienceandRemoteSensing, 2002, 40(12):2553-2566.

[10]ZhaXJ,FuRS,DaiZY,etal.Noisereductionininterferogramsusingthewaveletpackettransformandwienerfiltering[J].IEEEGeoscienceandRemoteSensingLetters, 2008, 5(3):404-408.

[11]BuadesA,CollB,MorelJM.Anon-localalgorithmforimagedenoising[C]//ProceedingsofIEEEInternationalConferenceonComputerVisionandPatternRecognition,SanDiego:IEEE, 2005.

[12]DabovK,FoiA,KatkovnikV,etal.Imagedenoisingbysparse3-Dtransform-domaincollaborativefiltering[J].IEEETransactionsonImageProcessing, 2007, 16(8): 2080-2095.

[13]MaggioniM,BoracchiG,FoiA,etal.Videodenoising,deblockingandenhancementthroughseparable4-Dnonlocalspatiotemporaltransforms[J].IEEETransactionsonImageProcessing, 2012, 21(9): 3952-3966.

[14]TasdizenT.Principalneighborhooddictionariesfornonlocalmeansimagedenoising[J].IEEETransactionsonImageProcessing, 2009, 18(12):2649-2660.

[15]CoupeP,YgerP,BarillotC,etal.Fastnonlocalmeansdenoisingfor3DMRimages[C]//ProceedingsofInternationalConferenceonMedicalImageComputingandComputerAssistedIntervention,Berlin:Springer, 2006.

[16]BoulangerJ,KervrannC,BouthemyP,etal.Patch-basednonlocalfunctionalfordenoisingfluorescencemicroscopyimagesequences[J].IEEETransactionsonMedicalImaging, 2010, 29(2): 442-454.

[17]ChatterjeeP,MilanfarP.Patch-basednear-optimalimagedenoising[J].IEEETransactionsonImageProcessing, 2012, 21(4):1635-1649.

[18]LeeJS,HoppelKW,MangoSA.IntensityandphasestatisticsofmultilookpolarimetricandinterferometricSARimagery[J].IEEETransactionsonGeoscienceandRemoteSensing, 1994, 32(5):1017-1027.

[19]ChatterjeeP,MilanfarP.Isdenoisingdead?[J].IEEETransactionsonImageProcessing, 2010, 19(4):895-911.

[20]ChatterjeeP,MilanfarP.Practicalboundsonimagedenoising:fromestimationtoinformation[J].IEEETransactionsonImageProcessing, 2011, 20(5): 1221-1233.

[21]BianY,BryanM.InterferometricSARphasefilteringinthewaveletdomainusingsimultaneousdetectionandestimation[J].IEEETransactionsonGeoscienceandRemoteSensing, 2011, 49(4):1396-1416.

[22]GhigliaDC,PrittMD.Two-dimensionalphaseunwrapping:theory,algorithmsandsoftware[M].USA:Wiley, 1998.

[23]HanCM,GuoHD,WangCL,etal.EdgepreservingfilterforSARimages[J].ChineseHighTechnologyLetters, 2003,7(13): 11-15.

[24]WangZ,BovikAC,SheikhHR.Imagequalityassessment:fromerrorvisibilitytostructuralsimilarity[J].IEEETransactionsonImageProcessing, 2004, 13(4):600-612.

Modified patch-based locally optimal wiener for interferometric phase filtering

WANG Yang1, HUANG Haifeng1, DONG Zhen1, WU Manqing1,2

Basedontheinterferometricphasefilteringproblemofsyntheticapertureradar,amodifiedpatch-basedlocallyoptimalwieneralgorithmwasproposed.TheproposedalgorithmwasthelinearminimummeansquareerrorestimatorundertheGaussianadditivenoiseconditionandjointlyestimatedthefirstmomentandsecondmomentoftheimage,namely,meanandcovarianceusingnon-localmeanswhichwasthestate-of-arttechnique.Whenappliedtointerferometricphasefiltering,twomodificationswereproposedaccordingtothespatialvariationofthenoise.First,meanvalue,insteadofmedianvalue,wasusedintheestimationofthenoisestandarddeviation.Second,thenumberofclusterswasdeterminedadaptivelyaccordingtotheratioofthemaximumvaluetothemeanvalueofthenoisestandarddeviation.Experimentalresultsonbothsimulationandrealdatashowthatthemodifiedpatch-basedlocallyoptimalwieneriseffectiveandissuperiortotheotherthreealgorithms.

interferometricphasefiltering;patch-basedlocallyoptimalwiener;linearminimummeansquareestimation

2014-10-11

國家自然科學基金資助項目(91438202)

汪洋(1983—)，男，四川綿陽人，博士研究生，E-mail：wycbx8384@163.com；黃海風(通信作者)，男，副研究員，博士，碩士生導師，E-mail：haifeng0728@vip.sina.com

10.11887/j.cn.201504017

http://journal.nudt.edu.cn

TN

A