改进分块局部最佳维纳滤波算法的干涉相位滤波*

汪 洋，黄海风，董 臻，吴曼青,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