王 爽 于佳平 刘 坤 侯 彪 焦李成

(西安电子科技大学智能感知与图像理解教育部重点实验室 西安 710071)

基于双边滤波的极化SAR相干斑抑制

王 爽*于佳平 刘 坤 侯 彪 焦李成

(西安电子科技大学智能感知与图像理解教育部重点实验室 西安 710071)

在对极化SAR数据进行相干斑抑制过程中，保持其极化信息是首要考虑的问题，而对极化SAR数据各元素进行独立的滤波是不可取。依据极化SAR数据形式及噪声模型，该文将双边滤波推广至极化SAR数据处理，推导出一种新的相似性距离度量公式，所提方法可以直接对极化协方差矩阵C或极化相干矩阵T进行处理。该方法能够保持极化信息，有效地抑制相干斑，同时能够更好地保持点目标、纹理结构等，在处理大幅面的极化 SAR数据时简单、快速、有效。

遥感图像；极化SAR(Pol-SAR)；双边滤波；距离度量；协方差矩阵C

1 引言

极化合成孔径雷达(Polarimetric Synthetic Aperture Radar, Pol-SAR)是遥感探测的重要手段，但是由于其相干成像原理，极化SAR图像中存在着严重的相干斑噪声[1]，影响了数据信息的提取以及后续的解译和分类。因此，斑点噪声的抑制是极化SAR数据处理中的重要问题，也是近几年国内外雷达图像处理研究的热点之一。

当图像中存在明显的边缘纹理信息时，通常的邻域运算会将边缘模糊化，对于高频细节的保护效果不明显。双边滤波是一种非线性的滤波方法，结合了图像的空间邻近度和像素相似性，同时考虑了空域信息和灰度信息，不仅能够有效地去除图像的噪声，还能很好地保持边缘纹理信息。

双边滤波由于其优越性得到了更广泛应用[2－5]。2009年，西安电子科技大学智能感知与图像理解教育部重点实验室Jiao等人[6]将双边滤波扩展到遥感SAR图像的相干斑抑制，2011年，Yang等人[7]对此进行了改进。2013年，Hondt等人[8]将双边滤波应用于极化SAR数据图像，该算法是基于黎曼距离度量的极化SAR相干斑抑制，其存在的问题是算法假设极化数据的统计分布是符合高斯分布，对数据进行了对数变换，导致结果会产生一定的辐射偏差。

本文依据极化SAR数据的噪声模型及复矩阵，推导出了一种新的相似性距离度量公式，将双边滤波应用于极化SAR数据的相干斑抑制。所提出方法可以直接对极化协方差矩阵C或极化相干矩阵T进行处理，能够很好地保持极化信息。

2 双边滤波的理论基础

1998年，Tomasi和Manduchi最先提出了双边滤波(Bilateral filtering)[9]的概念，它是在基于空间分布的高斯滤波函数基础上提出的。该方法结合了图像的空间邻近度和像素相似性，将二者非线性结合求得权重系数值，进而自适应滤波后获得平滑低噪图像。对于边缘附近的像素不会受到较远像素太多的影响，保证了对边缘附近像素值的保存。

普通自然图像的空域双边滤波算法[9]：

式中I(x)表示原像素值，表示滤波后的像素值，Ω表示邻域图像块的大小，fs(⋅)和fr(⋅)分别表示高斯核函数下的空间邻近距离及像素的相似距离，其表达分别为：

式中σs和σr是两个平滑参数，通常取经验值，K表示归一化函数，其公式为：

将双边滤波的表达式改写成权重系数的形式为：

其中，wi就是熟知的权重系数，其公式是：

双边滤波的权重系数是由Gauss权重系数和图像的亮度信息共同构成的。σr的加入可以对σs的变化做出相应的补偿，例如当σs相对来说比较大时，结合的像素就会相对比较多，导致输出的图像变得模糊，由于σr的限制，那些亮度与中心像素相差较大时，相应的权重较小，同时还去除了高频噪声。与其它局部算法相比较，其优点是不仅可以保证图像滤波效果，保持图像的边缘细节信息不失真，而且不需要进行复杂的运算，尤其是当窗口较大时，能够极大降低滤波过程的时间复杂性。

3 极化SAR图像的双边滤波

由单极化 SAR的噪声模型推广得到的极化SAR乘性噪声模型得到了广泛的应用[10－13]。对于单视极化SAR，雷达在非匀质区域所接收到的电磁波中包含着非匀质区域的空间波动所引起的纹理特性和相干斑噪声两个因素[14]。

极化域的乘性相干斑模型为：

式中，y表示雷达所接收到的回波矢量，v表示相干斑噪声矢量，t表示标量纹理变量。

对于多视极化SAR，由于

通常认为纹理变量的空间相关性比极化相干斑噪声特征矢量的空间相关性更高，因此可认为用于多视平均的像元的纹理变量是相等的，从而有

其中V代表相干斑的多视协方差矩阵，式(3)就是多视情况下极化乘性相干斑噪声模型，此时，Y服从一定条件的Wishart分布[15]。

目前大多数极化SAR数据的相干斑抑制、地物分类、目标识别都是在极化协方差矩阵C或极化相干矩阵T[16]的基础上进行处理的。

目标的极化协方差矩阵C定义为：

目标的极化相干矩阵T定义为：

则极化相干矩阵T可表示为：

极化协方差矩阵C和极化相干矩阵T可以相互转化，其转化公式为：

因此，将双边滤波应用于极化SAR图像的降斑中，需要处理的不再是单一的向量而是极化协方差矩阵C或极化相干矩阵T。我们将式(8)推导为基于极化矩阵权重系数的计算公式：

当处理的极化数据为极化协方差矩阵C时，Σ=C，同样，当处理的极化数据为极化相干矩阵T时，Σ=T。式中wi为权重系数，表达式为：

将双边滤波扩展到极化SAR数据上，其重点在于寻求一种有效的相似性距离度量方法，进一步来求解fr[⋅]得到相干斑抑制后的数据图像。

通过第3节的介绍可知，双边滤波最初提出是用于去除自然图像中的加性噪声，该方法采用了简单有效的欧式距离进行像素相似性度量：

式中，I(⋅)表示像素值。

对于乘性噪声模型的图像，采用比值距离可以度量像素之间的相似性

为了符合极化SAR的数据形式及噪声模型，满足当两个矩阵元素越相似，获得的权重越大，并且使其具有仿射不变性，针对上述问题，本文推导了一种新的相似性度量公式——改进的比值距离度量公式：

根据极化SAR数据的相干斑统计分布及模型，本文所提算法应用于极化 SAR的相干斑抑制主要从以下几点出发。首先，极化SAR数据的噪声模型是乘性模型，改进的比值距离度量公式能更好地度量像素之间的相似性。其次，采用双边滤波的思想，同时考虑了图像的空间邻近度和像素相似性，不仅能够有效地去除图像的噪声，还能很好地保持边缘纹理信息。再次，本文所构造的距离公式具有仿射不变性，对极化数据没有进行变换处理，不会产生辐射偏差。最后，本文推导的公式是直接对极化SAR的极化相干矩阵T或极化协方差矩阵C进行处理，能够保持其极化信息。

本文所提出的双边滤波算法本质上是一种非线性局部空域滤波算法，平滑参数的选取对滤波结果有较大的影响。当δr数值选取太大时会导致模糊，边缘细节信息丢失，数值太小，则导致区域滤波不平滑，去噪效果不明显。图1给出了δr选取不同数值，实测数据加拿大Ottawa地区迭代2次时的实验结果图，图2为不同平滑参数对应区域A的等效视数曲线图。

可以看出随着δr数值不断地增加，等效视数也会增加，但增加的幅度会降低，等效视数越大，表明对同质区域的滤波越平滑。根据实验结果，综合考虑对细节的保持和斑点的去除效果，所给出的δr的经验值为3.1。

其具体的实验步骤如下：

步骤1 读取极化SAR数据，参数设置，邻域窗大小w=7×7，两个平滑参数σs= 2.9和σr= 3.1，迭代次数IT(IT＜5)；

步骤2 求解空间距离度量的权重系数；

步骤3 取极化SAR数据的一个像素矩阵Σ(x)，以像素矩阵Σ(x)为中心取7×7大小的邻域窗，记为图像块X；

步骤4 求解邻域图像块X内的相似性散度距离度量的权重系数；

图1 不同平滑参数对应的结果图Fig. 1 The result figure with different smoothing parameters

步骤5 将极化SAR数据的每一像素矩阵对应的图像块，进行上述步骤4处理；

步骤6 得到初步滤波后结果；

步骤7 在基本的去噪结果上，对步骤3～步骤6进行迭代处理，直到满足要求，得到最终的相干斑抑制结果；

步骤8 用Pauli向量法对滤波后极化协方差矩阵C生成伪彩色图。

图2 不同平滑参数对应区域A的等效视数曲线图Fig. 2 Influence of smoothing parameters on the ENL value in A area

4 实验结果及分析

为了对本文提出的相似性距离度量的有效性进行分析，论文对4幅真实极化SAR数据进行了实验。实验所采用的4幅真实极化SAR数据中，3幅为多视极化数据，1幅为单视极化数据。第 1幅为美国San Francisco金门大桥区域，如图3(a)所示，该数据是2000年AIRSAR获取的，分辨率为10 m×10 m，4视，图中有明显的海洋、森林、城区3类主要地物；第2幅为加拿大Ottawa地区，Convair雷达数据，10视，大小为222×342，如图4(a)所示；第3幅为荷兰Flevoland区域， 1989年AIRSAR机载数据，分辨率为12.1 m×6.7 m，4视，截取其中的256×256的数据块，如图5(a)所示；第4幅图像为中国西安地区，2010年RadarSAT-2数据，分辨率为8 m×8 m，单视，截取其中的256×256的数据块，如图6(a)所示。

采用的对比算法为：经典的空域局部滤波算法精致极化Lee滤波[17]与IDAN(intensity-driven adaptiveneighborhood)[18]，Chen等人[11]所提的非局部均值滤波，以及 D’Hondt等人[14]使用的基于 Kullback-Leibler距离、Riemannian距离的极化SAR双边滤波。

降斑评价指标采用主观视觉判断和客观评价指标，客观评价标准主要用同质区域的等效视数(ENL)，无参考图像空域质量评价(BRISQUE)[19,20]。极化信息的保持是极化 SAR数据相干斑抑制重点关注的问题，实验中采用散射保持率与极化特征比值图来进行评价。

等效视数是极化 SAR图像降斑领域比较公认的降斑评价指标。等效视数表示滤波后极化SAR图像同质区域的平滑程度或者斑点噪声残留程度，其值越大，平滑程度越高，相应的斑点噪声的抑制越好。在计算等效视数时，要选择尽可能大的同质区域作为测试数据。极化SAR采用Pauli分解后的伪彩图计算ENL，可由同质区域的标准差与均值的比值获得。

2012年，Mittal等人[19]提出了无参考图像空域质量评价(BRISQUE)的概念，2013年，Torres等人[20]将该指标应用于极化SAR的相干斑抑制评价。BRISQUE的优势是无需参考比对图像，无需采用变换域，不需要计算图像的振铃、扭曲、模糊等特性，通过计算图像局部正规化亮度系数统计信息来评价图像潜在的失真度，从而在空域上评价整幅图像的质量。因此论文在这里也将其作为客观评价指标。BRISQUE的评价值在[0 100]区间内，数值越小说明图像的质量越好，根据文献[19]该指标的计算公式为：

图3-图5分别是San Francisco, Ottowa, levoland地区数据的相干斑抑制结果图，表1-表3是整幅图像的 BRISQUE数值及手动所选择的匀质区域的ENL数值，实验给出了本文算法迭代2次与迭代4次的结果。

图6是中国西安地区的相干斑抑制结果图。对单视数据，由于非局部方法中相似性阈值的问题，作者所提供的程序无法得到滤波结果，因此没有给出非局部均值滤波的结果，表 4是整幅图像的BRISQUE数值。

本文算法在不同的极化数据上迭代2次便能得到较好的斑点抑制效果，BRISQUE的值相对文中其它对比算法是最优的，对边缘细节纹理的保持和点目标的保持都也优于其它算法。但从同质区域的等效视数评价来说，迭代2次时该算法对匀质区域的滤波还存在一定的不足，算法迭代4次后，能够使同质区域更加平滑一致结果更理想，此时的边缘纹理的保持和点目标的保持仍然优于其它算法。

图3 San Francisco地区噪声抑制结果图Fig. 3 The speckle reduction results on the San Francisco image

图4 Ottowa地区噪声抑制结果图Fig. 4 The speckle reduction results on the Ottowa image

图5 Flevoland地区噪声抑制结果图Fig. 5 The speckle reduction results on the Flevoland image

表1 San Francisco地区的ENL与BRISQUE统计Tab. 1 ENL and BRISQUE values on the San Francisco image

表2 Ottowa地区的ENL与BRISQUE统计Tab. 2 ENL and BRISQUE values on the Ottowa image

表3 Flevoland地区的ENL与BRISQUE统计Tab. 3 ENL and BRISQUE values on the Flevoland image

图6 中国西安地区噪声抑制结果图Fig. 6 The speckle reduction results on the Xi’an image

表4 西安地区BRISQUE统计Tab. 4 BRISQUE values on the Xi’an image

由实验结果可以看出，在不同极化SAR数据中本文所提出的斑点抑制算法结果优于基于Kullback-Leibler距离和Riemannian距离的双边滤波算法。图5的数据由于影像中存在大量点呈现为极低散射回波点，造成部分低散射回波点作为地物细节信息被错误地保留下来。论文所提出的改进的比值距离方法和其它算法相比，虽然依然存在部分残留噪声引起的黑点，但是已经对上述问题有了较大的改善。

为了更好地分析本文方法对数据极化信息保持的效果，实验中采用散射保持率与极化特征比值图来进行评价。散射保持率的定义为滤波前后的表面散射(Ps)、体散射(Pv)、2次散射功率(Pd)值与原始功率值的比值：

从San Francisco极化SAR数据中挑选3种典型散射类型数据各100个，计算不同滤波后的散射保持率，最后对100个点的散射特性保持率求平均，数据在表5中给出。从表中可以发现，本文所提方法表面散射、体散射和2次散射都是保持最好的，其原因在于本文算法能更好地排除杂点。

表5 散射保持率Tab. 5 Scattering retention rate

从本文算法迭代2次后的San Francisco数据中选择海洋和森林区域中的100个点的平均极化特征信息做出极化特征图，并将不同滤波算法后的极化特征图与原数据的特征图做比值。图7是海洋区域的共极化特征比值图，图8是森林区域的交叉极化特征比值图。坐标Γ与坐标χ表示两个相位角度函数，取值范围分别为[0,90]和[0,180]。为滤波后极化信息保持得越好，极化特征图与原极化特征图越相同，理想的特征比值图像是在Z=1的平面上。从实验结果图可以看出，本文算法海洋区域共极化特征比值图所在的平面值为0.878最接近1，森林区域交叉极化特征比值图所在的平面值为0.882仍然最接近 1，并且其极化特征图像比值图出现异点幅度相对较少。因此，本文算法能够很好地保持其极化信息。

由以上实验结果可以看出，本文算法的滤波效果在一定程度上能与非局部均值相比，相比较非局部均值的时间复杂度，本文算法的更大优势在于减低了时间复杂度。非局部均值滤波算法的时间复杂度为。本文算法的时间复杂度为表6是对加拿大Ottawa地区的极化SAR数据进行算法运行时间对比，实验环境为MATLAB R2010a, Intel(R) Pentium(R) 1 CPU 2.4 GHz, Window XP Professional。

图7 San Francisco海洋区域共极化特征比值图Fig. 7 The ratio of co-polarization signatures in sea area of San Francisco

图8 San Francisco森林区域交叉极化特征比值图Fig. 8 The ratio of cross-polarization signatures in forest area of San Francisco

表6 Ottawa地区运行时间(s)Fig. 6 Running time of Ottawa area (s)

因此，本文算法能够达到对极化SAR数据图像降斑的要求，具有以下优势：(1)能够得到很好的相干斑抑制效果；(2)对整个极化协方差矩阵C矩阵或极化相干矩阵T直接进行处理，能够很好地保持极化信息。(3)同时考虑了图像的空间信息与亮度信息，能够更好地保持边缘纹理特征及点目标。(4)本文算法简单快速、易于理解。

5 结束语

本文根据极化SAR数据的特性，提出了一种新的相似性度量公式，从而将双边滤波应用于极化SAR数据相干斑抑制。该算法是直接对整个极化协方差矩阵C或极化相干矩阵T进行处理的。通过对真实极化SAR数据的实验仿真表明，该算法无论是对匀质区域的滤波还是边缘纹理、点目标和极化信息的保持上都达到对噪声抑制的要求，并且具有一定的优势，同时该算法简单、快速，可以在有效时间内处理大幅面的极化SAR数据。

[1] Kostinski, Alexander B, Boerner,et al.. On foundations of radar polarimetry[J].IEEE Transactions on AntennasPropagation, 1986, 34(12): 1395-1404.

[2] Rydell J, Knutsson H, and Borga M. Bilateral filtering of fMRI data[J].IEEE Journal of Selected Topics in Signal Processing, 2008, 2(6): 891-896.

[3] Zhang M and Gunturk B K. Multiresolution bilateral filtering for image denoising[J].IEEE Transactions on Image Processing, 2008, 17(12): 2324-2333.

[4] Anand C S and Sahambi J S. MRI denoising using bilateral filter in redundant wavelet domain[C]. IEEE Region 10 Conference, TENCON, Kyoto, Japan, 2008: 1-6.

[5] Paris S and Durand F. A Fast Approximation of the Bilateral Filter Using a Signal Processing Approach[M]. Springer Berlin Heidelberg: Computer Vision-ECCV 2006, 2006: 568-580.

[6] Zhang W G, Liu F, and Jiao L C. SAR image despeckling via bilateral filtering[J].Electronics Letters, 2009, 45(15): 781-783.

[7] Zhang W G, Zhang Q, and Yang C S. Improved bilateral filtering for SAR image despeckling[J].Electronics Letters, 2011, 47(4): 286-288.

[8] D’Hondt O, Guillaso S, and Hellwich O. Iterative bilateral filtering of polarimetric SAR data[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2013, 6(3): 1628-1639.

[9] Tomasi C and Manduchi R. Bilateral filtering for gray and color images[C]. IEEE Sixth International Conference on Computer Vision, Bombay, India, 1998: 839-846.

[10] Lee J S, Ainsworth T L, and Chen K S. Speckle filtering of dual-polarization and polarimetric SAR data based on improved sigma filter[C]. IEEE International Geoscience and Remote Sensing Symposium, Boston, MA, USA, 2008, 4: IV-21-IV-24.

[11] Chen J, Chen Y, An W,et al.. Nonlocal filtering for polarimetric sar data: a pretest approach[J].IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(5): 1744-1754.

[12] Lee J S, Grunes M R, Schuler D L,et al.. Scattering-modelbased speckle filtering of polarimetric SAR data[J].IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(1): 176-187.

[13] Liu G, Huang S, Torre A,et al.. The multilook polarimetric whitening filter (MPWF) for intensity speckle reduction in polarimetric SAR images[J].IEEE Transactions on Geoscience and Remote Sensing, 1998, 36(3): 1016-1020.

[14] 王超, 张红. 全极化合成孔径雷达图像处理[M]. 北京: 科学出版社, 2008: 60-61.

Wang Chao and Zhang Hong. Full Polarization Synthetic Aperture Radar Image Processing[M]. Beijing: Sciences Press, 2008: 60-61.

[15] Lopes A and Séry F. Optimal speckle reduction for the product model in multilook polarimetric SAR imagery and the Wishart distribution[J].IEEE Transactions on Geoscience and Remote Sensing, 1997, 35(3): 632-647.

[16] 王超, 张红. 全极化合成孔径雷达图像处理[M]. 北京: 科学出版社, 2008: 15-16. Wang Chao and Zhang Hong. Full Polarization Synthetic Aperture Radar Image Processing[M]. Beijing: Sciences Press, 2008: 15-16.

[17] Lee J S, Grunes M R, and de Grandi G. Polarimetric SAR speckle filtering and its implication for classification[J].IEEETransactions on Geoscience and Remote Sensing, 1999, 37(5): 2363-2373.

[18] 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.

[19] Mittal A, Moorthy A, and Bovik A. No-reference image quality assessment in the spatial domain[J].IEEE Transactions on Image Processing, 2012, 21(12): 4695-4708.

[20] Torres L, Sant’Anna S J S, da Costa Freitas C,et al.. Speckle reduction in polarimetric SAR imagery with stochastic distances and nonlocal means[J].Pattern Recognition, 2014, 47(1): 141-157.

[21] Yamaguchi Y, Moriyama T, Ishido M,et al.. Four-component scattering model for polarimetric SAR image decomposition[J].IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(8): 1699-1706.

王 爽(1978-)，女，2007年于西安电子科技大学获得工学博士学位；现为西安电子科技大学智能感知与图像理解教育部重点实验室教授，博士生导师，主要从事SAR/POLSAR处理与分析、稀疏表示、机器学习等方面的研究工作。

E-mail: shwang@mail.xidian.edu.cn

于佳平(1988-)，女，西安电子科技大学硕士研究生，研究方向为雷达图像处理。

刘 坤(1985-)，男，西安电子科技大学博士研究生，主要研究领域为 SAR图像处理、极化 SAR图像处理、机器学习等。

E-mail: lk314159@163.com

侯 彪(1974-)，男，西安电子科技大学教授，博士生导师，智能感知与图像理解教育部重点实验室副主任，IEEE会员，IET西安分会执行委员会委员，中国电子学会高级会员，陕西信号处理学会理事，教育部创新团队成员。主要研究方向为遥感图像解译、压缩感知、稀疏表示等。

E-mail: avcodec@163.com

焦李成(1959-)，男，西安电子科技大学教授，博士生导师，智能感知与图像理解教育部重点实验室主任，教育部创新团队负责人，主要从事智能信息处理、SAR影像处理与解译、压缩感知与稀疏表示等方面的研究工作。

E-mail: lchjiao@mail.xidian.edu.cn

Polarimetric SAR Speckle Reduction Based on Bilateral Filtering

Wang Shuang Yu Jia-ping Liu Kun Hou Biao Jiao Li-cheng
(Key Laboratory of Intelligent Perception and Image Understanding of Ministry of Education, Xidian University, Xi’an 710071, China)

The priority during speckle reduction in Polarimetric Synthetic Aperture Radar (Pol-SAR) data is to maintain the polarization information. For this reason, it is undesirable to separate each element of polarimetric SAR. Based on the multiplicative noise model and complex matrix, an effective distance similarity measure method is derived, which can process covariance matrixCor coherent matrixTdirectly. Experiments show that this method can effectively reduce speckle and maintain the polarization information, point targets, and texture structure. Moreover, it can handle large polarimetric SAR with simplicity and ease.

Remote sensing image; Polarimetric SAR (Pol-SAR); Bilateral filtering; Distance measure; Polarimetric covariance matrixC

中国分类号：TN957

A

2095-283X(2014)01-0035-10

10.3724/SP.J.1300.2014.13133

2013-12-18收到，2014-03-23改回；2014-03-27网络优先出版

61173092, 61271302)，新世纪优秀人才支持计划(NCET-11-0692)，陕西省科学技术研究发展计划项目 (B07048)，以及教育部“长江学者和创新团队发展计划”(IRT1170)资助课题*通信作者： 王爽 shwang@mail.xidian.edu.cn