基于极化SAR相干斑模型的非局部均值相干斑抑制

杨国辉，王 爽，刘立军，吕 龙，刘轶群

(1.兰州大学 信息科学与工程学院,甘肃 兰州 730013； 2.武警工程大学 密码工程学院,陕西 西安 710086； 3.西安电子科技大学 人工智能学院，陕西 西安 710071)

0 引言

极化合成孔径雷达(Polarimetric Synthetic Aperture Radar，Pol-SAR)是一种多参数、多通道的雷达成像系统[1-2]，是在单极化合成孔径雷达(Synthetic Aperture Radar，SAR)基础上发展而来。相比单极化SAR[3-5]，Pol-SAR根据接收和发射电磁波的方式不同，而具有不同的极化方式。通过测量目标每个分辨单元内的散射回波信息获得极化复散射矩阵S，目标复散射矩阵中包含了更多的散射和极化信息，为图像的理解、目标检测、目标识别和地物分类提供了更多的信息[6-7]。但是由于特殊的相干成像机理，极化SAR数据存在固有的相干斑噪声[8-9]，严重影响了后续特征提取、地物分类和目标检测等。因此，Pol-SAR相干斑噪声抑制依然是Pol-SAR数据处理与分析的关键问题，仍然是国内外合成孔径雷达图像处理与研究的重要方向之一。

早期的Pol-SAR数据相干斑抑制有两类，其中一类包括Novak等人[10]提出的极化白化滤波器,J.S.Lee[11]提出的最优加权滤波和极化矢量滤波[12]等，均是采取各通道间的极化相关信息通过数据融合处理达到相干斑抑制的效果；后期提出的算法基本上是上述算法的修正与改进，该类算法能够兼顾抑制相干斑和保持图像的分辨率，但是由于不同极化通道之间的串扰而不能很好地保持原始数据的极化散射特性，进而影响到后续的图像处理应用。另一类滤波方法主要借鉴了自然图像去噪和SAR图像相干斑抑制的思想，该类算法主要包括基于均值滤波的Boxcar滤波、基于Lee 滤波的Refined Lee滤波[11]以及Sigma滤波和改进Sigma滤波[13]等，但这类算法的不足是不能在有效抑制相干斑噪声的同时保持地物主散射特性和边缘纹理特征。针对上述算法的不足，如何在去除相干斑噪声的同时又能保留极化相关性和边缘特征等细节信息的问题，许多研究者们开展了边缘保持滤波算法的研究，提出的算法主要有：各向异性扩散法[14]、非局部均值滤波算法[15]和双边滤波算法[16]等。各向异性扩散法需要求解偏微分方程，但由于某些给定初始值上的病态性问题会导致处理过程不稳定，因此对高梯度值脉冲噪声的图像去噪效果不佳。相比而言，非局部均值滤波算法对于图像纹理和结构信息保持具有很好的效果，广泛应用于图像的去噪处理中，缺点是逐像素计算块相似度，计算量大，难以满足实时处理的需要。

本文详细介绍了Pol-SAR数据的相干斑噪声模型，提出了基于极化SAR相干斑模型的形状自适应非局部均值的相干斑去噪方法。本文方法可以直接对极化协方差矩阵C进行处理，根据对角线元素和非对角线元素噪声模型的不同分别进行处理，处理结果能够很好地保持极化相关信息。

1 极化SAR相干斑模型

1.1 极化SAR的数据描述

Pol-SAR数据由散射矩阵S来简单表示，散射矩阵S是一个2×2的复矩阵，包含了散射目标的相位和幅度信息[17]，具体如下：

(1)

式中，h，v分别表示水平极化和垂直极化；Shv表示水平极化接收和垂直极化发射，为其对应雷达照亮区域目标的散射系数，，其余元素表示类似。

对于Pol-SAR数据的处理与分析，为了便于表征，常采取将目标散射矩阵矢量化，在后向散射约定中，根据互易定理有Shv=Svh，因此目标向量K简化为：

(2)

假设Pol-SAR数据服从高斯分布，目标矢量k为一个服从多元零均值复高斯概率分布的随机变量，高斯概率密度函数(PDF)如下：

(3)

式中，H为复向量的共轭转置；C为协方差矩阵，表示如下：

[C]=E{KKH}=

(4)

单极化SAR数据的多视处理经常是平均邻近几个像素点以降低斑点噪声。在统计遍历性和局部同质性假设下，由空间平均替代统计估计得到协方差矩阵Z如下：

(5)

式中，L为视数，它是由平均像素斑点噪声来降低斑点噪声[18]。这种情况下，相干斑噪声的降低的同时降低了空间分辨率。多视协方差矩阵的Z统计PDF如下[20]：

(6)

根据极化SAR数据描述，极化SAR的每个分辨单元通过一个二维复散射矩阵S表示(见式(1))，将散射矩阵S向量化得到目标散射向量K，由目标向量的共轭积得到协方差矩阵，Pol-SAR数据间的关系如图1所示。

(7)

式中，z表示辐射幅度；φ表示辐射相位差。为了获得式(7)的相干斑模型，首先从复平面推导出相位差φ的噪声模型，然后将幅度信息z添加进去，最后得到完整的相干斑噪声模型。

图1 极化SAR数据关系

文献[19]研究表明，根据式(3)中数据分布的假设，极化SAR数据相干斑噪声模型包含乘性和加性噪声，该模型也被推广为多视极化SAR数据，其中，C对角线元素为乘性噪声，非对角线元素幅度部分为乘性噪声、相位部分为加性噪声。

1.2 极化SAR的相干斑模型

通常将极化SAR相干斑模型认为与单极化SAR相干斑模型相似，未考虑极化SAR数据特征，直接将极化SAR相干斑看成是乘性噪声。2003年，C.Lopez-Martine博士根据Pol-SAR的数据特征，系统地论述了Pol-SAR相干斑模型，并提出了乘-加性相干斑噪声模型[19]。本文基于该模型，设计Pol-SAR相干斑抑制算法。

S可以完全表征散射过程中确定的目标，但对于分布式的散射，它不能描述分布散射体即散射过程、随机指标等[19]，只能用C来描述，因此根据式(7)，Pol-SAR数据的相干斑模型可以通过C的进行推理，C的各元素的噪声模型表示如下[20]：

exp(jφx)+ψ(nar+jnai)，

(8)

式中，ψ表示在2个通道的平均能量：

(9)

|ρ|exp(jφx) 是复相关系数，表征2个极化通道Spq和Srs之间的相关性，表达式如下：

(10)

复相关系数的幅度与重要的滤波参数Nc 是一致的，Nc 表示为：

(11)

(12)

nm是乘性相干斑噪声的组成部分，其表征乘性噪声的均值和方差如下：

E{nm}=1，

(13)

var{nm}=1 。

(14)

nar+jnai是加性相干斑噪声，nar表示实部，nai表示虚部，其表征加性噪声模型的均值和方差如下：

E{nar}=E{nai}=0，

(15)

(16)

式(8)表示了2个通道相干斑噪声来源于乘性噪声源nm和加性噪声源na(na=nar+jnai)。由式(10)可以看出，复相关系数是确定乘性噪声和加性噪重要的参数，因此在相干斑抑制时，必须考虑复相关系数才能保证相关性的保持和更好的噪声抑制。

根据式(8)，单视和多视极化SAR数据的协方差矩阵噪声模型可简单表示为[20]：

Z=C+Nm+Na，

(17)

式中，矩阵Nm包含的是乘性噪声元素存在于整个协方差矩阵，而矩阵Na包含的是加性噪声信息，仅存在于C的非对角线元素相位上，因此在非对角线元素去噪时，要进行噪声分离，将幅度上的乘性噪声分离出来再处理。

2 形状自适应非局部均值滤波

2.1 非局部均值滤波

与传统空域滤波相比，非局部均值滤波[15](NL-means)是利用图像的冗余信息和某种相似性度量，在非局部窗口内或整幅图像内寻找像素块之间的相似性进行加权滤波。通过在整幅图像中寻找，从局部邻域的思想发展到了非局部，开创了图像去噪的新局面。 NL-means滤波数学表达式如下：

(18)

(19)

2.2 形状自适应非局部均值滤波

非局部均值滤波的相似窗是以像素为中心的方形窗，因此在边界处或纹理复杂的区域难以找到更多的相似块，因此导致信息利用不充分。本文提出了形状自适应非局部均值滤波算法，选取不同形状的相似窗有利于保持图像局部几何特性，使用各种掩膜形状相似窗来代替普通正方形相似块像素邻域，有利于信息的利用和保持。

本文提出了形状自适应非局部均值滤波算法，为了利用任意形状形似块的信息，其相似度距离公式表示如下：

(20)

式中，P(τ)表示所需的第τ类形状，如图2所示，本文用了8类不同形状的相似窗；x，y表示图像Z中相似块的中心像素。

图2 不同形状的形状掩膜Fig.2 Shape masks of different shapes

图2中每个形状掩膜中的黑色表示在该位置上的掩码为0，灰色表示在该位置上的掩码为1。

式(20)是基于加性噪声模型提出的相似度距离公式，对于极化SAR 的协方差矩阵而言，其对角线元素是乘性噪声模型，为了适用于乘性噪声模型，西安电子科技大学的钟桦等人提出了在贝叶斯框架下的乘性相干斑噪声的非局部均值相似度测量函数[21]：

(21)

根据式(20)和式(21)推导出适用于乘性噪声模型的形状自适应非局部均值滤波的相似度函数：

(22)

式中，P(k)表示第k类形状。

按照非局部均值滤波的式(19)，权值函数w(x,y)表示为：

(23)

上式是针对乘性噪声的非局部均值的权重函数，h为可调节滤波参数。

针对乘性噪声模型的形状自适应非局部均值的数学表达为：

(24)

图像中的任意像素x，可由K个不同形状相似邻域的估计得到：

(25)

由式(25)可以看出，对于像素x的估计值有N个滤波结果(本文N=8)。从N个结果中找出一个最优结果或者线性组合作为最终滤波结果，为了简化，本文对N个结果取其线性加权平均：

(26)

3 Pol-SAR相干斑抑制

图3 对角线元素的相干斑抑制原理Fig.3 Speckle filtering principle of diagonal elements

然后，按照非对角线元素的去噪原理，首先进行对角线元素去噪，然后计算复相关系数ρ:

(27)

式中，m，n为待处理图像尺寸大小。复相关系数是指通道Sp，Sq之间的相关性。

将得到复相关系数ρ的模值|ρ|和相位φ，对其相位exp(jφ)部分进行欧拉变换：

exp(jφ)=cos(φ)+jsin(φ)。

(28)

图4 非对角线元素去噪原理Fig.4 Denoising principle of off-diagonal elements

偏差因子B如下：

(29)

整个极化SAR协方差矩阵C的降斑的总体框架如图5所示。

图5 极化SAR整个协方差矩阵C降斑的总体框架Fig.5 Overall framework of despeckling of the entire covariance matrix C of polarimetric SAR

由图5可以看出，对于整个极化SAR协防矩阵C的元素要对角线和非对角线元素分别进行处理，尤其非对角线元素去噪处理更为复杂，也是极化SAR相干斑抑制的难点。

4 实验结果与评估

通过2组实验，验证了本文算法的有效性。将本文提出的滤波算法与其他经典滤波器相比，，采用评价指标EPD-ROA 和ENL对不同滤波算法性能进行客观评价，同时通过主观视觉效果进行评价。

4.1 客观评价指标

(1) 等效视数

等效视数[21](Equivalent Number of Looks，ENL)是衡量相干斑强度的参数，表示滤波后极化SAR图像上同质区域的平滑程度，其值越大说明平滑效果越好，也就是斑点噪声的抑制效果越好，定义如下：

(30)

式中，Imean，Iσ分别表示选定同质均匀区域的均值与方差；Hc是一个与极化SAR有关的常数，对于强度图像Hc=1，幅度图像Hc=4/π-1。本文中，首先在协方差矩阵C中进行相干斑的抑制处理，然后通过sinclair向量法将极化SAR数据生成span强度图像[22]，在span强度图上计算同质区域的等效视数来衡量同质均匀区域的相干斑抑制效果。

(2) 基于平均比率的边缘保持度量(EPD-ROA)

如何进行边缘的保持是极化SAR图像相干斑抑制需要考虑的一个重要因素，是度量降斑算法有效性的另一个重要指标。2011年，西安电子科技大学凤宏晓博士提出了基于平均比率的边缘保持度量[23](Edge-Preservation Degree based on of Average-EPD),定义如下：

(31)

式中，M是整幅图像像素总数；D1(i)，D2(i)分别表示沿某一方向降斑图像上的相邻像素；O1(i)，O2(i)分别表示沿某一方向原始图像上的相邻像素；X，Y分别表示去噪后的图像和原始图像。EPD-ROA 可以用来衡量图像滤波后的边缘保持度，其值越接近1，说明边缘保持度越好。本文在协方差矩阵C中进行相干斑的抑制处理，然后通过sinclair向量法将极化SAR数据生成span强度图像，使用span强度图像上计算EPD-ROA。

(3) 极化散射特性保持

极化散射特性的保持是极化SAR相干斑抑制评价的一个重要方面，也是一种相对主观的极化SAR 图像去噪质量评价指标，可从2个方面定性分析：① 通过分析滤波算法本身可以定性分析，滤波算法进行设计时是否考虑各通道之间的相关性，相关性的保持有利于极化散射特性的保持；② 极化SAR相干斑去噪是为了后续目标分类、目标检测等应用服务，如极化散射特性保持良好，将有利于后续的应用，可根据去斑后分类的结果来评价散射特性保持是否良好。

4.2 实验验证

为了验证本文算法的有效性，实验中用2幅Pol-SAR 图像进行测试：① 图6(a)是加拿大Ottawa地区图像，来源于机载雷达的L波段的极化SAR图像(大小：222 pixel×342 pixel，简称：Ottawa)；② 图7(a)是荷兰Flevoland区域图像，来源于机载雷达获得的四视极化SAR图像(大小：256 pixel×256 pixel，简称：Flevoland)。实验中，本文算法(BM-NLM)分别与精致Lee 滤波(Re-Lee)[11]、增强Sigma滤波(Im-Sigma)[13]进行比较,其中Re-Lee滤波和Im-Sigma滤波均采用7×7滑窗。通过等效视数ENL、比值图像以及边缘保持度量ROA_EPD三个方面的客观评价指标对比，同时分析了极化散射特性保持情况。

(a) 原图

(b) Re-Lee

(c) Im-Sigma

(d) BM-NLM

(a) 原图

(b) Re-Lee

(c) Im-Sigma

(d) BM-NLM

对比各个算法的去噪平滑效果，在图6(a)和图7(a)中选取了3个同质区域，计算所选区域的ENL，结果如表1所示。

表1 等效视数比较

实验结果中，在Ottawa、Flevoland图中使用对比算法及本文算法的降斑结果的比值图评价和边缘保持度ROA_EPD评价结果，如表2所示。

表2 各滤波结果的比值图和边缘保持度评价

由表1和表2可以看出，本文算法与对比算法相比，能够很好去除同质均匀区域的相干斑噪声。在视觉效果上，从图6中右下方、图7中部以及图7中右上角同质均匀区域的去噪效果上看出，本文算法相比其他2种滤波算法平滑很多。对比显示，Im-Sigma滤波优于Re-Lee滤波，但Im-Sigma滤波仍然在边缘处残留的噪声较为明显，因此本文算法在极化SAR相干斑噪声抑制上具有明显优势。在边缘细节信息的保持上，表2中客观评价指标边缘保持度EPD-RO表明，大多数情况下本文算法具有最优的效果，在视觉效果上，本文算法不仅在平滑程度上有很好的效果而且在边缘信息保持上也有不错的效果。从图6和图7可以看出，本文算法在点目标、线目标及边缘等的保持上相比对于算法具有明显的优势。

各滤波后wishart分类结果如图8所示。

(a) 原始数据分类

(b) Re-Lee滤波后分类

(c) Im-Sigma滤波后分类

(d) BM-NLM滤波后分类

在极化散射特性保持上：

(1) 算法本身的分析：本文算法基于极化SAR相干斑模型(乘-加噪声模型)进行设计，相干斑模型基于极化数据特征推理得到。在计算非局部均值滤波权重时，采用span强度图进行计算复相关系数ρ，保持了极化相关性；

(2) 滤波后分类效果：由图8可以看出，将滤波后的数据采用wishart进行分类，本文算法相比Re-Lee、Im-Sigma后分类的结果好很多，充分说明在滤波的同时极化散射特性保持的良好。

综上所述，本文算法在极化SAR相干斑抑制上优势明显，不论在客观评价指标还是视觉上都有不错的效果，同时保持了极化散射特性，对极化SAR后续处理具有重要的意义。

5 结束语

本文基于极化SAR相干斑模型，提出了形状自适应非局部均值的相干斑噪声抑制方法。本文提出的去噪方法在抑制相干斑噪声的同时可以更有效地保持极化特性和边缘结构信息。从极化SAR相干斑抑制效果可以看出，该方法在视觉效果和客观评价指标上都具有优越性。当前，针对极化SAR相干斑抑制问题，还缺少对其质量评价体系深入而系统的研究，相应的评价体系对于极化SAR相干斑抑制算法设计及极化SAR后续应用处理也至关重要，后续将进一步针对极化散射特性保持进行更深入的研究。