CN109035152A - 一种合成孔径雷达图像非局部均值滤波方法 - Google Patents

Abstract

本发明提出了一种合成孔径雷达图像非局部均值滤波方法，属于雷达图像处理领域。本发明将异性质测度引入非局部均值算法，改进了非局部均值算法的权值测量方式。首先，针对传统的非局部均值方法在计算距离测度采用等权值加权的方法，本文引入变差系数计算图像块窗口内像素点距离测度的加权权值，由于变差系数能够表征图像的纹理方向，因此引入变差系数加权后的图像块相似度能够更加有效的捕捉在纹理和方向上具有相似性的图像块，从而避免边缘性和方向性的模糊；然后，设计了一种基于变差系数的自适应退化参数函数，自适应的调整退化参数，可以有效的保证异质区得到保护不被过度平滑。所以本发明在滤除相干斑的同时，能够更有效的保持纹理边缘信息。

Description

一种合成孔径雷达图像非局部均值滤波方法

技术领域

本发明属于雷达图像处理领域，特别涉及一种合成孔径雷达图像非局部均值滤波方法。

背景技术

合成孔径雷达(Synthetic Aperture radar,SAR)是一种主动式的高分辨成像传感器，能够全天时、全天候的进行观测，穿透性强，广泛应用于包括军事侦查、灾情监测和地表覆盖层探测等领域。但是由于SAR成像是利用回波相干叠加进行成像，导致SAR图像不可避免的存在相干斑噪声，这些相干斑噪声降低了SAR图像的目视效果，对后续的图像分割和目标识别也带来困难。因此，SAR 图像的相干斑抑制对SAR图像的分割、分类和识别具有重要的意义。

SAR图像的相干斑滤波方法可以分为空域滤波和变换域滤波。空域滤波指将待滤波像素的邻域像素加权平均得到滤波后的灰度值，这类方法的代表有Lee 滤波、Frost滤波、Gamma-Map滤波等，他们都是基于特定的相干斑分布模型提出的，因此只能对特定分布的SAR图像有较好的滤波效果。变换域滤波指将图像从时域通过傅里叶变换或者小波变换等方法变换到相应的变换域，然后在变换域做滤波处理，最后逆变换回时域二维图像的方法，这类方法利用了多尺度、多方向的分析方法，常见的方法有小波变换、曲波变换、剪切波变换等。这两类方法从本质上来说，都是基于局部或者半局部的滤波方法，而没有考虑图像本身的结构冗余信息。

为了充分利用图像自身的结构冗余信息，提出了非局部均值(Non-Local Means,NLM)算法，该算法利用图像的结构冗余，在图像内寻找与目标图像块相似的图像块，并根据图像块间的相似性对像素点赋予相应的权值，利用加权平均的方法求得目标像素点的估计值，该方法在滤除噪声的同时可以很好的保持图像的纹理信息，因此也得到了越来越多的学者关注，并在NLM算法的基础上进行了改进。现有技术中，还有一方法，其在NLM算法的基础上，结合小波域滤波和维纳滤波的思想，提出了BM3D滤波方法，该方法分为基础估计和最终估计两次估计，初步估计把相似的块聚集成三维数组，对三维数组采用联合硬阈值进行滤波，把滤波后的图像块返回到他们原始的位置，得到重叠的块估计，然后加权平均得到初步估计的图像，最终估计是利用块相似性分别把含噪图像和初步估计的图像聚集成两个三维图像，然后使用联合维纳滤波，把滤波后的图像块放到原始位置得到重叠估计图像，加权后得到最终估计图像，取得了不错的效果。但是上述方法都只能适用于加性高斯白噪声，而不能在具有乘性噪声的SAR图像上，因此，很多学者对NLM算法在SAR图像上的应用进行了研究。因此，提出了基于图像块概率分布的(probability patch-based,PPB)非局部均值滤波方法，该方法利用了图像的噪声分布特性，在计算块相似性的时候采用了最大似然估计，并且在迭代过程中估计会越来越准确，该方法在SAR图像上取得了很好的滤波效果；其外，另一种方法证明了比率测度能够更稳健的表征SAR 图像的距离测度，并提出基于比率测度的匀质区非局部均值滤波方法，但是该方法只适用于匀质区的SAR图像去噪；在BM3D的基础上还提出了一种 SAR-BM3D方法，改进块相似性度量的方法，并解决了小波系数收缩的问题，在SAR图像上取得了非常好的滤波效果；还有一种方法，在比率测度非局部均值滤波的基础上，对距离测度的窗口权值做了改进，使Hongxiao Feng提出的比率测度非局部均值滤波能够推广到异质区。但是这些方法在计算距离测度时，使用的都是窗口内像素点灰度距离的均值，或者是高斯加权值，而没有考虑图像的纹理方向信息，而SAR图像中存在着大量的纹理方向信息，这些信息应该成为距离测度的一个重要依据。此外，在NLM算法中，权值是依赖于距离测度的高斯核函数，其中退化参数对滤波效果有非常重要的影响，然而很少有文献对退化参数进行研究，一般是根据经验值进行设置，没有实现自适应去噪的效果。

发明内容

为了解决上述问题，本发明提出了一种合成孔径雷达图像非局部均值滤波方法，引入变差系数(coefficient of variation,CV)，构造一种根据图像自适应变化的各向异性的二维加权窗；此外，提出一种基于变差系数的自适应调整退化参数的方法，使图像在同质区域有很好的去噪效果，同时在异质区能够保持边缘纹理细节。

一种合成孔径雷达图像非局部均值滤波方法，包括以下步骤：

步骤1，获取合成孔径雷达幅度图像，初始化参数，其中，初始化的参数包括搜查窗大小Ds、图像窗大小ds、退化参数系数γ及S曲线倾斜度参数ξ；

步骤2，计算得到所述图像的相干斑分布，基于所述相干斑分布得到相干斑的噪声方差；

步骤3，根据所述图像的标准差和均值计算图像变差系数；

步骤4，根据所述图像变差系数，计算自适应退化参数；

步骤5，基于所述变差系数，构造自适应滤波窗口；

步骤6，基于异性质测度的非局部算法，得到像素点的值。

进一步地，所述步骤1还包括以下流程：

设定像素点i为位于图像位置i的像素，图像块i为以像素点i为中心的大小为ds的图像块，搜索区Ω_i为以像素点i为中心的大小为Ds的图像块。

进一步地，所述步骤2包括以下流程：

当图像像素间距和雷达成像分辨单元相当时，噪声功率是非相关的，相干斑是非相干的乘性噪声，表示为

υ＝u·n

其中，υ为观测值，u为雷达反射系数(RCS)，n为相干斑噪声，根据图像的标准差和均值可得图像的等效视数

其中，L为等效视数，σ为图像的标准差，μ为图像的均值，相干斑完全发展，其强度条件分布函数为

其中，Γ(·)为Gamma函数，设定雷达反射系数为1，得到完全发展的相干斑噪声分布

计算得到噪声强度的方差为

所述图像的幅度值为强度的平方根，可得观测值的幅度条件分布为

设定雷达反射系数为1，可得幅度噪声的分布函数

计算得到噪声方差为

进一步地，所述步骤3包括以下流程：

变差系数为

其中，CV为像素点的局部变差系数。

进一步地，所述步骤4包括以下流程：

基于S曲线以自适应调整变差系数，S曲线表示为

其中，参数τ用于控制S曲线的偏移量，参数ξ用于控制S曲线的倾斜度；

根据所述变差系数，对于图像中的每个像素点，自适应退化参数表示为

其中，h′为自适应退化因子，h＝γσ，CV_max为整个图像的最大局部变差系数。

进一步地，所述步骤5包括以下流程：

根据所述步骤2中的乘性噪声模型，且噪声是均值为单位值的高斯白噪声，一个图像块的变差系数表示为

其中，Ψ为观测图像υ中的一个图像块，若该图像块在同质区，则 E(u_Ψ ²)＝[E(u_Ψ)]²，变差系数可简化为

根据噪声值为单位值，可得在同质区的图像块变差系数近似为图像噪声方差，根据图像局部变差系数构造的各向异性的滤波窗口表示为

A_CV(i,k)＝exp(-|CV(i,k)-CV(i)|²)

其中，A_CV(i,k)为以i为中心的各向异性的加权窗，CV(i,k)为以i为中心的图像块中第k个像素点的变差系数。

进一步地，所述步骤6包括以下流程：

步骤61，计算图像块相似性度量

基于比率测度和所述自适应滤波窗口，相似性测度表示为

其中，d_CV(i,j)为以像素点i和j为中心的图像块的相似性测度，Ψ为图像块，Z_CV(i)为图像块中变差系数的归一化因子，函数Ξ(·)的表达式为Ξ(x)＝x+1/x， u_i(k)为以i为中心的图像块的第k个像素点的灰度值，u_j(k)为以j为中心的图像块的第k个像素点的灰度值；

步骤62，计算像素点权值

根据NLM算法，基于高斯核函数，计算像素点间的权值

其中，为归一化系数；

步骤63，计算像素点像素值

根据NLM算法，对于像素点i，将其搜索区Ω_i中的像素点按权值叠加，得到像素点i的估计值

其中，u(j)为像素j对应的灰度值。

本发明的有益效果：本发明将异性质测度引入非局部均值算法，改进了非局部均值算法的权值测量方式。首先，针对传统的非局部均值方法在计算距离测度采用等权值加权的方法，本文引入变差系数计算图像块窗口内像素点距离测度的加权权值，由于变差系数能够表征图像的纹理方向，因此引入变差系数加权后的图像块相似度能够更加有效的捕捉在纹理和方向上具有相似性的图像块，从而避免边缘性和方向性的模糊；然后，设计了一种基于变差系数的自适应退化参数函数，自适应的调整退化参数，可以有效的保证异质区得到保护不被过度平滑。所以本发明在滤除相干斑的同时，能够更有效的保持纹理边缘信息。

附图说明

图1为本发明的流程图。

图2为图1中步骤6的流程图。

具体实施方式

下面结合附图对本发明的实施例做进一步的说明。

请参阅图1，本发明提出的一种合成孔径雷达图像非局部均值滤波方法，通过以下步骤实现：

步骤1，获取合成孔径雷达幅度图像，初始化参数，其中，初始化的参数包括搜查窗大小Ds、图像窗大小ds、退化参数系数γ及S曲线倾斜度参数ξ。

本实施例中，输入的图像SAR机载幅度图，图像中包括同质区图像块A 和异质区图像块B，图像块大小为30×30，用于图像滤波结果评估，波段为X 波段，分辨率为3米，图像大小为256×256。初始化的参数包括搜查窗大小Ds、图像窗大小ds、退化参数系数γ及S曲线倾斜度参数ξ。

本实施例中，设定像素点i为位于图像位置i的像素，图像块i为以像素点 i为中心的大小为ds的图像块，搜索区Ω_i为以像素点i为中心的大小为Ds的图像块。

步骤2，计算得到所述图像的相干斑分布，基于所述相干斑分布得到相干斑的噪声方差。

本实施例中，步骤2通过以下流程实现：

当图像像素间距和雷达成像分辨单元相当时，噪声功率是非相关的，可以认为相干斑是非相干的乘性噪声，表示为

υ＝u·n

其中，υ为观测值，u为真实的雷达反射系数(RCS)，n为相干斑噪声，根据输入图像的标准差和均值可得图像的等效视数，其表达式为

其中，L为等效视数，σ为图像的标准差，μ为图像的均值，对于视数为 L的SAR图像，若相干斑完全发展，其强度条件分布函数可以表示为

其中，Γ(·)为Gamma函数，设定雷达反射系数为1，得到完全发展的相干斑噪声分布

计算得到噪声强度的方差为

SAR图像的幅度值是强度的平方根，可得观测值的幅度条件分布为

设定雷达反射系数为1，可得幅度噪声的分布函数

计算得到噪声方差为

步骤3，根据所述图像的标准差和均值计算图像变差系数。

本实施例中，步骤3通过以下流程实现：

变差系数是图像标准差和均值的比值，即对于SAR图像中的每一个像素，该像素位置的局部变差系数就是取该像素点的局部邻域图像块来估计该像素的变差系数，像素点的局部邻域大小和图像块大小定义一致，取以像素点为中心的ds×ds的图像块用于计算中心像素点的变差系数：

其中，σ_Ψ和μ_Ψ为像素点局部邻域的标准差和均值。

步骤4，根据所述图像变差系数，计算自适应退化参数。

本实施例中，步骤4通过以下流程实现：

在传统的NLM算法中，根据相似性测度计算权值的指数函数的退化因子在整幅图中是恒定不变的，因此退化因子的选择十分重要，如果太大，则会导致过度平滑，不利于图像的纹理边缘的保持，如果太小，则在同质区不能达到满意的滤波效果，由于变差系数能够区分同质区和异质区，本发明结合变差系数，提出一种自适应的退化因子计算方法使退化因子能够自动调整，适应不同区域图像的滤波要求。由变差系数的定义可知，同质区的变差系数较小，接近于噪声方差，异质区的变差系数远大于噪声方差，结合上述对退化因子的分析，本发明利用S型曲线来自适应的调整变差系数，S型曲线表达式如下：

其中，参数τ用于控制S曲线的偏移量，参数ξ用于控制S曲线的倾斜度。

根据变差系数，对于图像中的每个像素点，自适应退化参数表示为

其中，h′为自适应退化因子，h和传统的非局部均值一样，参数h的取值和图像的噪声方差有关，是一个经验值，一般取h＝γσ，其中，γ取值范围在 [90,100]内可以取得较好的效果。由于同质区的变差系数接近噪声方差，异质区的变差系数远大于噪声方差，因此参数τ选择S型曲线的对称点位置等于噪声方差时计算得到的值，即在CV＝σ时，τ+ξ(CV/CV_max)＝0求解得到τ＝-ξ(CV/CV_max)，参数ξ用于控制S型曲线的倾斜度，经过大量实验，得到参数ξ的范围在[20,50]可以取得较好效果，要说明的是，这个步骤是针对整个图像的全局变差系数计算的，其中参数表示图像的全局最大变差系数。CV_max为整个图像的最大局部变差系数，用于对变差系数进行归一化。

步骤5，基于所述变差系数，构造自适应滤波窗口。

本实施例中，步骤5通过以下流程实现：

根据所述步骤2可知SAR图像的相干斑噪声可以被简化为乘性噪声，若噪声是均值为单位值的高斯白噪声，根据步骤B的乘性噪声模型，一个图像块的变差系数可以表示为

其中，Ψ为观测图像υ中的一个图像块，若该图像块在同质区，则 E(u_Ψ ²)＝[E(u_Ψ)]²，变差系数可简化为

因为噪声均值为单位值，因此可得也就是说因此在同质区的图像块变差系数可以近似为图像噪声方差；如果该图像块在异质区，那么 E(u_Ψ ²)＞＞[E(u_Ψ)]²，因此在异质区的图像块变差系数远大于图像噪声方差，也就是说，变差系数在边缘纹理处的值远大于同质区域，利用图像局部变差系数的这个特点，可以构造自适应的各向异性的二维滤波窗口，在距离测度中引入各向异性的滤波窗口有利于增大相似性测度的精度，可以更加有效的选择相似性高的像素点，从而能够更加有效地保护边缘纹理目标和点目标；同时，在同质区，因为变差系数近似于噪声方差，滤波窗口近似平坦，因此可以保证各向同性的滤波，也就是说，在不影响同质区域的情况下，充分利用了边缘纹理方向信息进行加权，有效的捕捉了图像的方向细节。根据图像局部变差系数构造的各向异性的滤波窗口可以表示为：

A_CV(i,k)＝exp(-|CV(i,k)-CV(i)|²)

其中，A_CV(i,k)为以i为中心的各向异性的加权窗，CV(i,k)为以i为中心的图像块中第k个像素点的变差系数。如果i位于纹理上，利用变差系数可以将边缘像素赋与高权值，而与边缘方向无关的像素点赋与低权值，从而达到方向选择的目的，那么A_CV(i,k)沿图像纹理有形成滤波窗口；如果i位于点目标上，点目标的变差系数远大于周围像素点，因此可以得到中心加权的形式，那么A_CV(i,k)分布类似于高斯滤波窗口；如果i位于同质区，由于同质区的变差系数都接近为噪声方差，因此各个像素点几乎赋与相同权值，那么A_CV(i,k)是一个分布均匀的平坦滤波窗口。

步骤6，基于异性质测度的非局部算法，得到像素点的值。

请参阅图2，步骤6通过以下流程实现：

步骤61，计算图像块相似性度量

根据相干斑可以简化建模为乘性噪声这个特点，本发明使用比率测度来计算相似性度量，结合步骤5构造的自适应滤波窗口，对于搜索区中以i和j为中心的图像块，相似性测度表示为

其中，d_CV(i,j)为以像素点i和j为中心的图像块的相似性测度，Ψ为图像块，为图像块中变差系数的归一化因子，函数Ξ(·)的表达式为Ξ(x)＝x+1/x，r(i,j)＝u_i(k)/u_j(k)为以i和j为中心的图像块的第k个像素点的灰度值u_i(k)和u_j(k)之间的比值，可以看出，r(i,j)越接近1，说明两个灰度值越接近，r(i,j)越远离1，说明两个灰度值相差越大，而函数Ξ(·)是一个在[0,+∞]先递减再递增的函数，并在r(i,j)等于1时取得最小值。

步骤62，计算像素点权值

根据NLM算法，使用高斯核函数计算像素点的权值，其中采用上述的图像块相似性测度作为距离，退化系数采用步骤4中计算的自适应退化系数，最后进行归一化即可，像素点间的权值可以表示为

其中，为归一化系数，h′为自适应退化因子，用于控制指数函数的衰减速度。

步骤63，计算像素点像素值

根据NLM算法，对于像素点i，将其搜索区Ω_i中的像素点按权值叠加，得到像素点i的估计值

其中，u(j)为像素j对应的灰度值。

本发明对每个像素点的像素值估计为NL_CV[u](i)。本方法滤除相干斑的同时，能够更加有效的保持纹理边缘信息。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (7)

1.一种合成孔径雷达图像非局部均值滤波方法，其特征在于，包括以下步骤：

步骤1，获取合成孔径雷达幅度图像，初始化参数，其中，初始化的参数包括搜查窗大小Ds、图像窗大小ds、退化参数系数γ及S曲线倾斜度参数ξ；

步骤2，计算得到所述图像的相干斑分布，基于所述相干斑分布得到相干斑的噪声方差；

步骤3，根据所述图像的标准差和均值计算图像变差系数；

步骤4，根据所述图像变差系数，计算自适应退化参数；

步骤5，基于所述变差系数，构造自适应滤波窗口；

步骤6，基于异性质测度的非局部算法，得到像素点的值。

2.如权利要求1所述的合成孔径雷达图像非局部均值滤波方法，其特征在于，所述步骤1还包括以下流程：

设定像素点i为位于图像位置i的像素，图像块i为以像素点i为中心的大小为ds的图像块，搜索区Ω_i为以像素点i为中心的大小为Ds的图像块。

3.如权利要求2所述的合成孔径雷达图像非局部均值滤波方法，其特征在于，所述步骤2包括以下流程：

当图像像素间距和雷达成像分辨单元相当时，噪声功率是非相关的，相干斑是非相干的乘性噪声，表示为

υ＝u·n

其中，υ为观测值，u为雷达反射系数(RCS)，n为相干斑噪声，根据图像的标准差和均值可得图像的等效视数

其中，L为等效视数，σ为图像的标准差，μ为图像的均值，相干斑完全发展，其强度条件分布函数为

其中，Γ(·)为Gamma函数，设定雷达反射系数为1，得到完全发展的相干斑噪声分布

计算得到噪声强度的方差为

所述图像的幅度值为强度的平方根，可得观测值的幅度条件分布为

设定雷达反射系数为1，可得幅度噪声的分布函数

计算得到噪声方差为

4.如权利要求3所述的合成孔径雷达图像非局部均值滤波方法，其特征在于，所述步骤3包括以下流程：

变差系数为

其中，CV为像素点的局部变差系数。

5.如权利要求4所述的合成孔径雷达图像非局部均值滤波方法，其特征在于，所述步骤4包括以下流程：

基于S曲线以自适应调整变差系数，S曲线表示为

其中，参数τ用于控制S曲线的偏移量，参数ξ用于控制S曲线的倾斜度；

根据所述变差系数，对于图像中的每个像素点，自适应退化参数表示为

其中，h′为自适应退化因子，h＝γσ，CV_max为整个图像的最大局部变差系数。

6.如权利要求5所述的合成孔径雷达图像非局部均值滤波方法，其特征在于，所述步骤5包括以下流程：

根据所述步骤2中的乘性噪声模型，且噪声是均值为单位值的高斯白噪声，一个图像块的变差系数表示为

其中，Ψ为观测图像υ中的一个图像块，若该图像块在同质区，则E(u_Ψ ²)＝[E(u_Ψ)]²，变差系数可简化为

根据噪声值为单位值，可得在同质区的图像块变差系数近似为图像噪声方差，根据图像局部变差系数构造的各向异性的滤波窗口表示为

A_CV(i,k)＝exp(-|CV(i,k)-CV(i)|²)

其中，A_CV(i,k)为以i为中心的各向异性的加权窗，CV(i,k)为以i为中心的图像块中第k个像素点的变差系数。

7.如权利要求6所述的合成孔径雷达图像非局部均值滤波方法，其特征在于，所述步骤6包括以下流程：

步骤61，计算图像块相似性度量

基于比率测度和所述自适应滤波窗口，相似性测度表示为

其中，d_CV(i,j)为以像素点i和j为中心的图像块的相似性测度，Ψ为图像块，Z_CV(i)为图像块中变差系数的归一化因子，函数Ξ(·)的表达式为Ξ(x)＝x+1/x，u_i(k)为以i为中心的图像块的第k个像素点的灰度值，u_j(k)为以j为中心的图像块的第k个像素点的灰度值；

步骤62，计算像素点权值

根据NLM算法，基于高斯核函数，计算像素点间的权值

其中，为归一化系数；

步骤63，计算像素点像素值

根据NLM算法，对于像素点i，将其搜索区Ω_i中的像素点按权值叠加，得到像素点i的估计值

其中，u(j)为像素j对应的灰度值。