基于非局部均值滤波的SAR图像抑噪

2020-07-20张文哲潘军蒋立军仲伟敬张雪峰范博文卞宇涛

世界地质 2020年2期

张文哲，潘军，蒋立军，仲伟敬，张雪峰，范博文, 卞宇涛

1.吉林大学地球探测科学与技术学院，长春 130026；2.中国人民解放军63771部队，西安 714000

0 引言

合成孔径雷达是一种可全天时、全天候持续观测的主动式微波遥感器，其在灾害评估、目标分类、目标检测识别和军事侦察等民用及军事领域均有重大的应用。然而，对于SAR图像而言其分辨尺寸是信号波长的几十倍或更多倍，因此每一分辨单元内所收到的信号是由多个散射点矢量叠加[1]。由于不同散射点到传感器距离不同，因此，同一地面单元内的散射点回波相位均不相同，这使得散射回波矢量叠加后获得的幅度图像的像元值并不完全取决于地物后向散射系数而是偏离正常幅度值随机起伏变化。由于这个原因，使得SAR图像上同质区域内产生与正常地物亮度差异明显且呈颗粒状随机分布的相干斑噪声,这些噪声的存在使人工目视判读及利用像元值来计算并获取信息的过程产生错误计算而获得错误结果，进而严重影响了后续的SAR图像的解译与应用过程[2]。因此，研究如何在有效地抑制相干斑噪声的同时保留更多纹理信息对于SAR图像分割、分类和自动目标识别等方面具有重要意义[3]。

SAR图像滤波的主要目的在于通过对图像噪声的抑制使均质区域更加平滑，使图像纹理区保留更多信息，对SAR图像的滤波方法研究较多的为空域滤波与变换域滤波两类[1]。空域方法利用了邻域像元之间的相关性，对固定大小窗口内的像元进行加权平均来得到中心像元的像元值，此种方法能够较好地去除相干斑噪声，但会使图像变得模糊而丢失纹理细节信息，常见方法有Lee滤波、Kuan滤波及Gamma MAP滤波等[1]。变换域类方法通过空间域与频率域间的变换来更好地抑制噪声，通过多方向分析来保留纹理结构，能够在抑噪的同时保留纹理细节。但此类方法涉及多个中间过程的计算，增加了算法的复杂性，且不同方法的适用性有限，此类方法以小波变换为主，有小波变换、平稳小波变换等方法。小波变换虽能更好保留纹理，但其计算更为复杂[4]，相比而言空域滤波显得更易于实现。

在空域滤波中，为在滤波的同时保留更多纹理结构信息，Buades et al.[5]在2005年提出针对加性白噪声的非局部均值(non-local means,NL-means)滤波算法。与以往算法利用单个像素点的相似性构造权重不同，该算法利用图像块间的相似性信息，在整个图像内寻找与目标图像块相似的图像块，依据图像块的相似性计算权重，加权平均求得像素的估计值[6]。NL-means算法在图像均质区域可以取得良好的抑噪效果，但因相似块的过加权而使纹理区过于平滑，未能有效地保留图像纹理信息。因此，有不少国内外学者开始在NL-means算法的基础上展开研究，并且取得了一定的成果。Deledalle et al.[7]提出在图像块概率分布(probabilistic patch-based，PPB)基础上进行加权平均的抑噪算法，通过真实SAR图像实验，此方法取得了较好的抑噪效果。刘书君等人[8]利用局部统计信息将图像划分为同质与异质区，并且使用不同方法计算两种区域内相似块间的权值，提高了相似块匹配精度，增强了抑噪性能。易子鳞等人[9]提出引入结构相似性指数来改进图像块间的相似性度量，此方法在抑制图像噪声的同时能够更好地保持图像纹理结构信息。朱磊等人[10]利用局部统计信息，以相似框与中心区域局部均值比代替高斯加权欧氏距离，以变差系数代替常系数，以此形成新的权值进行计算，较于传统算法具有更好的抑噪及边缘保持性能。以上算法存在处理过程繁琐或抑噪、纹理结构保留性能不佳的问题。

笔者从提高抑噪算法对图像纹理结构保留能力出发，采用两级滤波的方法，通过一级高斯滤波(Gauss filter)削弱图像噪声，以此来提高相似块度量及梯度计算的准确性。二级滤波采用改进的NL-means方法，将梯度引入相似性度量函数，保留图像纹理结构信息，结合二者提出了基于梯度的二级NL-means算法(gradient-based secondary non-local means，GS-NLM)。实验结果表明，所提算法在有效抑制相干斑噪声的同时能更好地保留图像纹理结构。

1 非局部均值算法(NL-means)

非局部均值算法[5]考虑到整幅图像中存在相似结构的信息，通过在图像中寻找与目标图像块具有相似结构的图像块来构建相似权重，最后通过加权平均来计算目标像元的滤波像元值。

算法可表述为：

NL[I](i)=∑j∈Ωw(i,j)I(j)

(1)

式中：NL[I](i)为滤波后的图像像元；I(j)为以像元j为中心的图像块中的像元；Ω为搜索窗内像元坐标域；w(i,j)权值大小由像元i与像元j的相似性决定，0≤w(i,j)≤1且满足∑jw(i,j)=1。

像元i与像元j的相似性取决于其灰度值向量I(Ni)与I(Nj)间的相似性，其中，Ni表示以像元坐标i为中心的方形邻域，二者相似性通过高斯加权的欧氏距离d(i,j)来计算，其表达式为：

(2)

式中：α表示高斯核标准差，且α>0。

欧氏距离越小则邻域灰度值向量越相似，相应的像元在加权平均中获得的权值越大，权重定义如下：

(3)

2 基于梯度的二级NL-means算法

NL-means算法在图像均质区抑噪效果明显，但在图像纹理区显得过为平滑，这是由于算法在相似块的权值计算上存在过加权问题，即只考虑了像元值的相似性而未考虑图像本身的结构信息。像元梯度在图像纹理区域变化明显，能够反映出图像纹理结构；在均质平坦区域梯度变化缓慢，可以描述图像的均匀连续性。因此，本文将图像梯度加入NL-means算法中图像块的相似性权值计算，提出基于梯度的二级NL-means(gradient-based secondary non-local means,GS-NLM)算法，通过综合计算相似块中对应像元间的像元值与梯度的相似性来提高算法在纹理边缘处的权值，以此提高算法对图像纹理结构信息的保持能力。

因NL-means算法是在图像噪声为加性高斯白噪声的基础上提出的[5]，而SAR图像噪声符合乘性模型。因此，首先需对SAR强度或幅度图像进行对数变换，将乘性相干斑噪声转换为加性高斯白噪声，然后对图像进行一级高斯滤波(Gauss filter)削弱图像噪声以此来减小噪声对下一步梯度计算的影响;其次用一级滤波的结果求取图像梯度,并用此梯度图像来改进权值计算对一级滤波图像进行二级改进的NL-means算法滤波;最后对图像进行指数变换以恢复SAR图像的原始辐射特性。本文GS-NLM算法流程如图1所示。

图1 GS-NLM算法流程图Fig.1 Flow chart of GS-NLM algorithm

2.1 一级高斯滤波

SAR图像的乘性噪声服从Gamma分布，对图像进行对数变换后，噪声转换为加性噪声，此时噪声分布可近似为高斯分布[11]。高斯滤波(Gauss filter)是一种线性平滑抑噪算法，它可以使图像平滑，减少图像因噪声存在而产生的明暗的锐利变化，该算法对于抑制服从正态分布的噪声效果较好。据滤波窗口内的像元与中心像元的空间欧氏距离，利用高斯函数来构建窗口内对应邻域像元的权值，则高斯滤波可表示为：

(4)

式中：Igf(x,y)为坐标(x,y)处的高斯滤波像元值；Ω为以坐标(x,y)处的像元为中心的m×m大小的方形窗口邻域；w(i,j)为坐标(x,y)处的像元权值。

2.2 二级改进NL-means算法

NL-means算法在图像块相似性度量中只考虑了图像块间的像元灰度相似性而未考虑图像本身的纹理结构信息，导致抑噪效果不佳，虽然能够较好地平滑图像但未能有效保留其纹理结构信息。

对于SAR图像，其纹理区域亮度变化明显，像元值变化较大，因此，加入对此变化反映明显的梯度，用像元灰度欧氏距离及梯度欧氏距离来对图像块进行相似性度量，使算法对图像纹理区更加敏感，滤波后的图像能够在平滑的同时保留其纹理结构信息。算法定义如下：

I′(x,y)=∑j∈Ωw(i,j)wg(i,j)Igf(j)

(5)

式中：I′(x,y)为坐标(x,y)处的滤波像元值；w(i,j)与Ω的含义同上文；

2.3 算法中一些参数取值的讨论

为减小图像噪声在图像块相似性度量和梯度计算时的负面影响，先对原始图像进行了一级高斯预滤波。为抑制噪声使图像不至于过度平滑而保留更多信息，在高斯滤波时选择大小为5×5像素的滤波窗口，经多次实验对比后将滤波参数值设置为0.6。

在式(3)中计算图像块间的相似性权值时，参数h控制指数函数的衰减速度，h越大指数函数衰减越慢，滤波图像越平滑，反之则指数函数衰减越快，滤波图像越达不到好的抑噪效果。在Buades et al.[5]的算法中，h≈10σ，即h的大小与图像噪声标准差成正比关系。对于本文实验所用的SAR幅度图像，采用Coupe et al.[12]提出的方法，即由图像噪声水平来确定平滑参数h(h≈2kσ2)，其中，k为一个常数，σ2为噪声方差。对于实验图像的噪声方差估计，采用Xie et al.[13]对SAR幅度图像噪声方差的估计方法。经本文对真实SAR图像的多次试验表明，当h≈0.38σ2时，抑噪效果最佳。对于控制梯度相似度指数函数衰减的参数，本文经多次试验取值为0.05。

3 实验结果与分析

3.1 实验设计

为验证本文所提出的算法的性能，对从欧洲航天局获得的ALOS PALSAR数据进行预处理,得到两幅单视真实SAR幅度图像，以此为基础进行试验。图像大小均为512×512像素，两幅图像分别记为SAR1和SAR2(图2)。

为对比显示本文算法在结构保持上的优越性，选取了经典的空域滤波算法：Refined Lee滤波算法、Gamma MAP滤波算法、双边滤波(BF)算法及NL-means算法来做比较。其中，增强Lee算法、Gamma MAP算法及双边滤波算法的滤波窗口大小均取7×7像素。为充分利用图像冗余信息的同时减小算法运算量，NL-means算法及本文算法的搜索窗口大小均选取21×21像素，相似图像块的窗口大小选取7×7像素[5]。

a.SAR1;b.SAR2。图2 实验所用真实SAR图像Fig.2 Real SAR images used in experiment

3.2 算法评价指标

对于滤波后的SAR图像抑噪质量的评价包括主观评价与客观评价。主观评价主要通过人眼对抑噪前后图像的对比观察来实现，客观评价指标有等效视数及结构相似度指数。本文同时采用主观评价与客观评价。

(1)等效视数

等效视数(equivalent number of looks，ENL)用于衡量相干斑噪声的抑制程度，表示为：

ENL=μ2/σ2

(6)

式中：μ表示图像像元灰度平均值；σ2表示对应的方差。ENL要在图像的均匀区域计算，本文选择在图像中手动选择一块区域，每种方法的抑噪结果均在此区域计算ENL的值，选取的均匀区域如图2白色框所示。ENL值越大，表明相干斑噪声程度越弱，抑噪效果越好。对于本文实验所用的幅度数据，ENL需乘以变差系数的平方，即4/π-1。

(2)结构相似度指数

结构相似度指数(structural similarity index，SSIM)[14]反应图像滤波前后的相似性，综合考虑了图像亮度、对比度和结构信息，其函数表示为：

(7)

式中：μi与σi、μj与σj分别代表图像i与图像j像元值的均值与方差；σij代表两图像间的协方差；C1与C2为两常数，通常取C1=(k1*L)2、C2=(k2*L)2；k1与k2为调节系数，且k1,k2≪1；L为图像灰度级。

(3)主观评价

主观评价包括对图像内点目标、线目标和纹理结构的保持程度，对图像的平滑效果及图像噪声的均匀与否等。

3.3 实验结果

将各个相干斑噪声抑制算法分别作用于两幅SAR图像，滤波结果如图3与图5所示，对应算法的比值噪声图如图4与图6所示。

图3与图5为不同算法对两幅SAR图像的抑噪结果，对比两图中的b-f可以看出：就抑噪程度看，Refined Lee、BF、NL-means、GS-NLM算法抑噪结果中斑点噪声较少，性能要优于Gamma MAP算法，而四者中NL-means与GS-NLM算法斑点噪声最少，故而抑噪性能最好；就图像平滑效果看，Gamma MAP算法远不够平滑，而Refined Lee与BF算法整体过于平滑导致图像模糊，NL-means与GS-NLM算法在图像均质区足够平滑且图像并不模糊；就图像纹理、点、线目标保持效果看，对比图像相同位置处的点线纹理，Gamma MAP算法保持最为显著，而GS-NLM算法次之且优于NL-means算法。

a.SAR1;b.Refined Lee;c.Gamma MAP;d.BF;e.NL-means;f.GS-NLM。图3 SAR1抑噪后的图像Fig.3 SAR1 images after despeckling

a. Refined Lee比值噪声图;b.Gamma MAP比值噪声图;c.BF比值噪声图;d.NL-means比值噪声图;e.GS-NLM比值噪声图。图4 SAR1抑噪后的比值噪声图Fig.4 SAR1 ratio noise images after despeckling

a.SAR1;b.Refined Lee;c.Gamma MAP;d.BF;e.NL-means;f.GS-NLM。图5 SAR2抑噪后的图像Fig.5 SAR2 images after despeckling

a. Refined Lee比值噪声图;b.Gamma MAP比值噪声图;c.BF比值噪声图;d.NL-means比值噪声图;e.GS-NLM比值噪声图。图6 SAR2抑噪后的比值噪声图Fig.6 SAR2 ratio noise images after despeckling

图4与图6为不同算法所计算出的比值噪声图，对比两图中的a-e可以看出，Gamma MAP算法噪声较为均匀，但其点、线处较光滑，说明纹理处未能有效去除噪声，在纹理处抑噪效果较差；Refined Lee与BF算法中含有较多图像点、线纹理且斑点噪声较少，说明抑噪能力不强，也未能保留多少图像纹理；NL-means与GS-NLM算法含有较少纹理且具有较多斑点噪声，说明抑噪能力较强且可更大程度地保留图像纹理。

表1与表2为不同算法抑噪结果的客观评价，对比可以看出：本文所提GS-NLM算法ENL为107.734 3与479.141 2，分别比其他算法高出3.7～72.5与18.5～346.2，而Gamma MAP算法值最小，说明本文算法对图像均质区抑噪平滑效果最好，而Gamma MAP算法最差；对比SSIM指数可以看出除Gamma MAP算法外，本文GS-NLM算法为0.792 8与0.750 1，均比其他算法大，说明除Gamma MAP算法外，本文算法滤波后的结果与原图像更相似，且纹理结构信息保留更多。