CN106780361A - 一种用于城市区域合成孔径雷达图像的噪声抑制方法 - Google Patents

Abstract

本发明公开了一种用于城市区域合成孔径雷达图像的噪声抑制方法，该方法包括如下步骤：对SAR图像实部/虚部数据上滑动窗口内样本数据估计Alpha‑Stable分布参数，记录所有滑动窗口的参数；对SAR图像幅度数据对数变换，在对数变换后图像上划分各对应于实部/虚部滑动窗口，计算各对应窗口内的Alpha‑Stable分布表达式和噪声分布表达式；计算Alpha‑Stable分布表达式与噪声分布表达式乘积的后验概率，取动态范围内最大后验概率的值作为滑动窗口中心象素的值；将所有经最大后验估计后的滑动窗口中心象素值作指数变换，得到滤波后幅度影像。本发明所述技术方案更有效地抑制城市区域SAR图像的斑点噪声，提高了SAR图像的质量。

Description

一种用于城市区域合成孔径雷达图像的噪声抑制方法

技术领域

本发明涉及一种合成孔径雷达图像噪声抑制方法，特别地涉及一种用于城市区域合成孔径雷达图像的噪声抑制方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar，简称SAR)全天候、全天时工作的特点，使合成孔径雷达遥感在固体地球科学、生态科学、水文科学和海洋科学等领域正发挥着越来越重要的作用，特别是近年来随着城市化进程的加快，合成孔径雷达凭借其成像优势在遥感应用中占有越来越重要的地位。但SAR图像特有的相干斑特别是城市区域相干斑噪声一直困扰着SAR图像的处理和解译等。

在过去的几十年中提出了许多SAR图像斑点噪声抑制方法，大致可分为两类：成像前的多视处理技术和成像后的处理技术。但多视处理为了提高辐射分辨率，忽略了不同的信息，如雷达散射截面随观测角的变化或者空间变化的影响，且图像的分辨率都会降低。成像后的处理技术SAR图像斑点噪声抑制方法，以空间域滤波技术为基础的相干斑抑制算法为主流，基于最大后验概率(Maximum a posteriori，简称MAP)，噪声在一定程度上被抑制，图像边缘的信息很少丢失，且无伪边缘产生，但此方法不适合于不同性质SAR图像区域(均匀区、非均匀区和孤立点)。

近年来采用更精确的Alpha-Stable分布模型代替传统的分布模型模拟SAR图像分布，以最大后验估计准则，应用于SAR图像斑点噪声领域取得显著效果。

发明内容

为了克服现有技术中的不足，本发明提供一种用于城市区域合成孔径雷达图像的噪声抑制方法，采用拟合程度更高的Alpha-Stable分布模型代替传统的统计分布模型，在处理城区SAR图像时，与几种常用的SAR图像滤波算法相比，在相干斑抑制和细节信息保持等方面具有明显的优势。

为了达到上述发明目的，解决其技术问题所采用的技术方案如下：

一种用于城市区域合成孔径雷达图像的噪声抑制方法，包括以下步骤：

步骤S1：对SAR图像实部/虚部数据上滑动窗口内的样本数据估计Alpha-Stable分布参数，记录所有滑动窗口的参数α_i,j、γ_i,j，其中，参数α和γ的下标i、j表示滑动窗口的行、列数；

步骤S2：对SAR图像幅度数据进行对数变换，在对数变换后的图像上划分各对应于实部/虚部的滑动窗口，计算各对应滑动窗口内的Alpha-Stable分布表达式：

其中，随机变量X服从阿尔法稳定分布的概率密度函数，X为SAR图像实部/虚部幅度值的自然对数，u∈[0,∞)，参数α∈(0,2]，决定该分布脉冲特性的程度，参数γ为分散系数，J₀(x)为第一类0阶贝赛尔函数；

和噪声分布表达式：

其中，随机变量N服从Nakagami分布的概率密度函数，N为SAR图像幅度值对数变换后加性斑点噪声，L是SAR图像的视数，Γ(F)是一个均值为1、方差为的具有Γ分布的的概率密度函数，F为SAR图像幅度值的乘性噪声幅度；

步骤S3：根据贝叶斯公式推出

并由最大后验概率准则：

即在X动态范围内使P_X(X)*P_N(N)值最大的X_i作为滑动窗口中心象素的值；

步骤S4：将所有经最大后验估计后的滑动窗口中心象素值，作指数变换得到滤波后幅度影像。

进一步的，所述步骤S1进一步包括如下子步骤：

步骤S1.1：采用滑动窗口，在图像各局部选取样本，按R²＝σ²/u²计算每一局部样本的相对标准差，其中，R、σ、u分别为SAR图像幅度值局部样本的相对标准差、标准方差、均值，构成相对标准差的直方图，并假设图像的相对标准差服从χ²分布；

步骤S1.2：根据所得的图像相对标准差直方图，估计χ²分布的参数，计算在所有的自由度f下的相对标准差的值R(f)；

步骤S1.3：采用滤波窗口邻域划分进行窗口中围绕中心像素的最大均匀区的搜索，即计算每个邻域的相对标准差Ci，剔除Ci>R(f)的邻域，余下的部分就是窗口内最大的均匀区，如果没有一个邻域使Ci≤R(f)，则缩小窗口尺寸，重新进行邻域划分和最大均匀区的搜索；

步骤S1.4：在搜索到的最大均匀区中采用均值滤波进行斑点抑制，其中搜索最大均匀区域的过程是以局部相对标准差为线索，通过变化窗口大小或采用不断地剔除窗口中非均匀的邻域得到最后的结果。

进一步的，步骤S1.1中所述图像各局部选取样本数在15～25个像素之间。

优选的，步骤S1.1中所述图像各局部选取样本数为8，所采用的滑动窗口为9×9。

进一步的，步骤S1.3中进一步包括计算滑动窗口内同质区域的均值u、标准差σ和相对标准差C_x：

当C_x≤C_u时，滤波区域为均匀区域，进行均值滤波；

当时，滤波区域为斑点完全生成区，进行Alpha-Stable分布最大后验概率滤波；

否则，剔除窗口内Ci最大的邻域，并假设窗口中剩下的所有邻域为同质区域。

本发明所述技术方案解决了SAR图像中抑制斑点噪声的问题，与传统的滤波方法相比，有以下两个关键特点:

一是采用一种改进的斑点噪声抑制算法，该算法以局部相对标准差为主要依据，将待处理像素归属为均匀区、斑点噪声完全生成区和孤立点目标，分别采用均值滤波、Alpha-Stable分布最大后验概率滤波和取中心象素值进行对应的处理，避免计算所有滑动窗口的Alpha-Stable分布，减少了计算量。二是引入边缘保持最优化(Edge PreservingOptimized Speckle，简称EPOS)滤波算法的思想，基于改变窗口形状的边缘保持最优化滤波算法，通过逐步剔除与中心像素不同质的邻域，然后应用最后剩余的邻域内点的统计量估算最后的结果，以达到更好地抑制斑点噪声，并保持图像的细节信息效果。本发明所述技术方案更有效地抑制城市区域SAR图像的斑点噪声，提高了SAR图像的质量。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单的介绍。显而易见，下面描述中的附图仅仅是本发明的一些实施例，对于本领域技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。附图中：

图1为本发明中基于Alpha-Stable分布的抑制斑点噪声流程图；

图2为本发明中基于Alpha-Stable分布可变窗口滤波算法流程图；

图3为本发明中EPOS算法窗口邻域划分示意图；

图4试验区域为选取的1993年11月13日2：32，ERS-1获取的上海市静安区苏州河附近256×256大小的原始幅度影像；(a)原始幅度影像(b)Kuan滤波(c)增强Lee滤波(d)本文滤波算法；

图5试验区域为2005年4月28日机载城区原始幅度影像，大小为256×256象素；(a)原始幅度影像(b)Kuan滤波(c)增强Lee滤波(d)本文滤波算法。

具体实施方式

以下将结合本发明的附图，对本发明实施例中的技术方案进行清楚、完整的描述和讨论，显然，这里所描述的仅仅是本发明的一部分实例，并不是全部的实例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。

如图1所示，本发明公开了一种用于城市区域合成孔径雷达图像的噪声抑制方法，包括以下步骤：

步骤S1：对SAR图像实部/虚部数据上滑动窗口内的样本数据估计Alpha-Stable分布参数，记录所有滑动窗口的参数α_i,j、γ_i,j，其中，参数α和γ的下标i、j表示滑动窗口的行、列数；

步骤S2：对SAR图像幅度数据进行对数变换，在对数变换后的图像上划分各对应于实部/虚部的滑动窗口，计算各对应滑动窗口内的Alpha-Stable分布表达式：

其中，随机变量X服从阿尔法稳定分布的概率密度函数，X为SAR图像实部/虚部幅度值的自然对数，u∈[0,∞)，参数α∈(0,2]，决定该分布脉冲特性的程度，参数γ为分散系数，J₀(x)为第一类0阶贝赛尔函数；

和噪声分布表达式：

其中，随机变量N服从Nakagami分布的概率密度函数，N为SAR图像幅度值对数变换后加性斑点噪声，L是SAR图像的视数，Γ(F)是一个均值为1、方差为的具有Γ分布的的概率密度函数，F为SAR图像幅度值的乘性噪声幅度；

步骤S3：根据贝叶斯公式推出

并由最大后验概率准则：

即在X动态范围内使P_X(X)*P_N(N)值最大的X_i作为滑动窗口中心象素的值；

步骤S4：将所有经最大后验估计后的滑动窗口中心象素值，作指数变换得到滤波后幅度影像。

进一步的，所述步骤S1进一步包括如下子步骤：

步骤S1.1：采用滑动窗口，在图像各局部选取样本，按R²＝σ²/u²计算每一局部样本的相对标准差，其中，R、σ、u分别为SAR图像幅度值局部样本的相对标准差、标准方差、均值，构成相对标准差的直方图，并假设图像的相对标准差服从χ²分布(Chi-squaredistribution，卡方分布)；

步骤S1.2：根据所得的图像相对标准差直方图，估计χ²分布的参数，计算在所有的自由度f下的相对标准差的值R(f)；

步骤S1.3：采用滤波窗口邻域划分进行窗口中围绕中心像素的最大均匀区的搜索，即计算每个邻域的相对标准差Ci，剔除Ci>R(f)的邻域，余下的部分就是窗口内最大的均匀区，如果没有一个邻域使Ci≤R(f)，则缩小窗口尺寸，重新进行邻域划分和最大均匀区的搜索；

步骤S1.4：在搜索到的最大均匀区中采用均值滤波进行斑点抑制，其中搜索最大均匀区域的过程是以局部相对标准差为线索，通过变化窗口大小或采用不断地剔除窗口中非均匀的邻域(变换滤波窗口形状)得到最后的结果。

进一步的，步骤S1.1中所述图像各局部选取样本数在15～25个像素之间。优选的，步骤S1.1中所述图像各局部选取样本数为8，所采用的滑动窗口为9×9。

进一步的，步骤S1.3中进一步包括计算滑动窗口内同质区域的均值u、标准差σ和相对标准差C_x：

当C_x≤C_u时，滤波区域为均匀区域，进行均值滤波；

当时，滤波区域为斑点完全生成区，进行Alpha-Stable分布最大后验概率滤波；

否则，剔除窗口内Ci最大的邻域，并假设窗口中剩下的所有邻域为同质区域。

本实施例中提供的用于合成孔径雷达图像的图像斑点噪声抑制算法，是利用VC++实现的，采用星载单视数SAR图像、机载SAR城区为试验区比较各滤波器在上的滤波效果。

具体实现步骤如下：

①对移动窗口(M×M)进行如图2所示的邻域划分，分别计算窗口内8个邻域相对于中心像素灰度值的相对标准差Ci(i＝1…8)，假设整个窗口属同质区域；

②计算窗口内同质区域的均值u、标准差σ和相对标准差C_x；

③如在C_x≤C_u时，滤波区域为均匀区域，进行均值滤波；如有可认为滤波区域为斑点完全生成区，进行Alpha-Stable分布最大后验概率滤波；否则，剔除窗口内Ci最大的那个邻域，并假设窗口中剩下的所有邻域为同质区域，返回步骤②，这里L为SAR图像的视数；

④在步骤③中进行Alpha-Stable分布最大后验概率滤波时，当估计的参数α＞2或者斑点完全生成区样本个数小于100个时，选择Kuan滤波。

⑤如窗口内邻域均被剔除，则缩小窗口尺寸(M＝M－2)并返回步骤①；

⑥如窗口尺寸已缩小到3×3，则认为该中心像素为点状目标，保留中心象素值。

算法的核心问题是在滤波滑动窗口中寻找与中心像素统计特性同属一类的像素点分布区域(滤波区域)，并选用适用算法在滤波区域进行滤波操作。该算法可实现滑动窗口尺寸和参加滤波像素的自动选择，在进行斑点噪声抑制的同时较大程度地保护了图像的边缘和细节信息。

ERS-1(欧洲遥感卫星1号卫星，特指欧洲遥感卫星1号卫星SAR载荷)上海城区各算法滤波前后的均值、均值偏移、标准差、等效视数统计结果分别如表1所示。

均值偏移度＝(滤波后图像均值－原始图像均值)/原始图像均值

表1上海城区SAR数据的滤波实验结果比较

	均值	均值偏移	标准差	等效视数ENL
					原始幅度影像	127.7	--	78.0	2.68
Kuan滤波	127.6	-0.078％	47.3	7.28
					增强Lee滤波	126.8	0.078％	37.1	11.67
本文滤波	125.8	-1.574％	31.3	17.76

城区机载SAR数据各算法滤波前后的均值、均值偏移、标准差、等效视数统计结果分别如表2所示。

表2机载城区SAR数据的滤波实验结果比较

	均值	均值偏移	标准差	等效视数ENL
					原始幅度影像	371.5	--	514.5	0.52
Kuan滤波	371.4	-0.027％	351.1	1.12
					增强Lee滤波	365.3	-1.669％	361.3	1.02
本文滤波	357.3	-3.773％	352.5	1.92

由图4、图5和表1、表2我们可以看到，Kuan滤波算法和增强Lee算法目视结果中存在一些影响视觉效果的亮点，本发明提出的基于Alpha-Stable分布可变窗口滤波算法对于ERS-1上海城区和机载城区SAR图像的相干斑滤除效果和边缘细节保持效果都较好，等效视数优于其它算法，但均值偏移、标准差指标上不如Kuan滤波算法和增强Lee算法。另外，基于Alpha-Stable分布可变窗口滤波算法，通过将子窗口划分为均匀区域、完全发育相干斑区域、强散射体区域，减少了需要计算Alpha-Stable分布的子窗口的个数，在计算效率上有一定的提高。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (5)

1.一种用于城市区域合成孔径雷达图像的噪声抑制方法，其特征在于，包括以下步骤：

步骤S1：对SAR图像实部/虚部数据上滑动窗口内的样本数据估计Alpha-Stable分布参数，记录所有滑动窗口的参数α_i,j、γ_i,j，其中，参数α和γ的下标i、j表示滑动窗口的行、列数；

步骤S2：对SAR图像幅度数据进行对数变换，在对数变换后的图像上划分各对应于实部/虚部的滑动窗口，计算各对应滑动窗口内的Alpha-Stable分布表达式：

$P_A\left(X\right)=e^{2X}{\∫}_0^{\mathrm{\∞}}u\exp (-{\mathrm{\γu}}^{\mathrm{\α}})J_0\left({\mathrm{ue}}^X\right)du$

其中，随机变量X服从阿尔法稳定分布的概率密度函数，X为SAR图像实部/虚部幅度值的自然对数，u∈[0,∞)，参数α∈(0,2]，决定该分布脉冲特性的程度，参数γ为分散系数，J₀(x)为第一类0阶贝赛尔函数；

和噪声分布表达式：

$p_A\left(N\right)=\frac{2L^L\exp (2NL-{\mathrm{Le}}^{2N})}{\mathrm{\Γ}\left(F\right)}$

其中，随机变量N服从Nakagami分布的概率密度函数，N为SAR图像幅度值对数变换后加性斑点噪声，L是SAR图像的视数，Γ(F)是一个均值为1、方差为的具有Γ分布的的概率密度函数，F为SAR图像幅度值的乘性噪声幅度；

步骤S3：根据贝叶斯公式推出

$P_{X|Y}\left(X\right|Y)=\frac{P_{Y|X}\left(Y\right|X)P_X\left(X\right)}{P_Y\left(Y\right)}=P_N\left(N\right)P_X\left(X\right)$

并由最大后验概率准则：

$\hat{X}\left(Y\right)=\arg {\mathrm{maxP}}_N\left(N\right)P_X\left(X\right)$

即在X动态范围内使P_X(X)*P_N(N)值最大的X_i作为滑动窗口中心象素的值；

步骤S4：将所有经最大后验估计后的滑动窗口中心象素值，作指数变换得到滤波后幅度影像。

2.根据权利要求1所述的一种用于城市区域合成孔径雷达图像的噪声抑制方法，其特征在于，所述步骤S1进一步包括如下子步骤：

步骤S1.1：采用滑动窗口，在图像各局部选取样本，按R²＝σ²/u²计算每一局部样本的相对标准差，其中，R、σ、u分别为SAR图像幅度值局部样本的相对标准差、标准方差、均值，构成相对标准差的直方图，并假设图像的相对标准差服从χ²分布；

步骤S1.2：根据所得的图像相对标准差直方图，估计χ²分布的参数，计算在所有的自由度f下的相对标准差的值R(f)；

步骤S1.3：采用滤波窗口邻域划分进行窗口中围绕中心像素的最大均匀区的搜索，即计算每个邻域的相对标准差Ci，剔除Ci>R(f)的邻域，余下的部分就是窗口内最大的均匀区，如果没有一个邻域使Ci≤R(f)，则缩小窗口尺寸，重新进行邻域划分和最大均匀区的搜索；

步骤S1.4：在搜索到的最大均匀区中采用均值滤波进行斑点抑制，其中搜索最大均匀区域的过程是以局部相对标准差为线索，通过变化窗口大小或采用不断地剔除窗口中非均匀的邻域得到最后的结果。

3.根据权利要求2所述的一种用于城市区域合成孔径雷达图像的噪声抑制方法，其特征在于，步骤S1.1中所述图像各局部选取样本数在15～25个像素之间。

4.根据权利要求3所述的一种用于城市区域合成孔径雷达图像的噪声抑制方法，其特征在于，步骤S1.1中所述图像各局部选取样本数为8，所采用的滑动窗口为9×9。

5.根据权利要求2所述的一种用于城市区域合成孔径雷达图像的噪声抑制方法，其特征在于，步骤S1.3中进一步包括计算滑动窗口内同质区域的均值u、标准差σ和相对标准差C_x：

当C_x≤C_u时，滤波区域为均匀区域，进行均值滤波；

当时，滤波区域为斑点完全生成区，进行Alpha-Stable分布最大后验概率滤波；

否则，剔除窗口内Ci最大的邻域，并假设窗口中剩下的所有邻域为同质区域。