CN103077508B - 基于变换域非局部和最小均方误差的sar图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于变换域非局部和最小均方误差的SAR图像去噪方法，主要解决现有技术对SAR图像去噪时产生边缘过平滑和点目标难以保持的问题。其实现步骤是：输入一幅SAR图像Y，并对其进行一层非下采样拉普拉斯分解，得到一幅低频图像Y_L和一幅高频图像Y_H；对Y_L利用PPB滤波器进行滤波，得到滤波后的图像对Y_H进行剪切波方向滤波器分解，得到各个方向的高频子带图像对利用高斯混合模型建模和MMSE估计进行去噪，得到去噪后的高频子带图对和进行逆剪切波变换，得到空域图像Y_Z；对Y_Z进行分类处理，得到最终去噪结果。本发明既能去除同质区域的噪声，也能很好地保持图像的边缘清晰，可用于对图像的预处理。

Description

基于变换域非局部和最小均方误差的SAR图像去噪方法

技术领域

本发明属于图像处理技术领域，涉及合成孔径雷达SAR图像的去噪，具体地说是一种基于变换域非局部和最小均方误差MMSE估计的SAR图像去噪方法，可用于图像的预处理阶段。

背景技术

合成孔径雷达具有全天时、全天候成像、高空分辨率和强穿透能力等特点，因而SAR图像在军事和民用方面均获得了广泛的应用。然而，由于系统的有限分辨率和相干性，SAR成像过程中总是不可避免地要产生一种斑点噪声。斑点噪声的存在严重影响了SAR图像的质量。SAR图像去噪需要对同质区域有效地抑制噪声，同时，还需要尽可能地保持边缘和点目标等重要信息不丢失，它是SAR图像处理和分析的关键步骤，也是保证后续图像处理如分割、目标识别等或图像解译等工作顺利完成的基础。

目前的SAR图像去噪算法主要在空间域和频率域进行，传统的基于空域的滤波算法有Lee滤波、Kuan滤波、Frost滤波和Gamma-MAP滤波等，这些基于空域的滤波方法，虽然都有相对较好的降斑结果，但难以保持图像中的边缘和细节信息。其原因在于它们都假设真实的SAR图像场景是平稳的，而这个假设在SAR图像的边缘和细节区域是不成立的。非局部NL均值滤波方法在自然图像去噪中取得了很好的结果，近些年在图像降噪等领域取得了广泛的应用。由于其简单的思想和良好的性能，在近几年也被推广到SAR图像降斑中。2009年Zhong H.等在文献“Bayesian Nonlocal means filterfor SAR image despeckling.2nd AsianPacific Conference on Synthetic ApertureRadar(APSAR),2009:1096-1099”中提出了一种基于Bayesian框架下的非局部均值对SAR图像降斑的算法，取得了不错的效果，同年C.A.Deledalle等在文献“IterativeWeighted Maximum Likelihood Denoising with Probabilistic Patch-Based Weights.IEEETransactions on Image Processing,2009,Vol.18,No.12,Pages:2661-2672”中提出了一种基于图像块的迭代非局部滤波器（Probabilistic Patch-Based PPB），它是非局部均值NL滤波器的一种扩展，该方法获得了非常好的降斑效果，被认为是目前性能最好的SAR图像降斑方法之一，但该方法仍会一定程度的模糊图像边缘。

基于变换域的去噪方法一般在小波域和后小波域进行，由于小波变换具有良好的多分辨分析和局部时频表示能力，对图像可以实现较好的稀疏表示，因而被广泛用于SAR图像去噪领域中。1997年Chipman H.A.等在文献“Adaptive bayesian waveletshrinkage.Journal of the American Statistical Association,1997,92(440):1413-1421”中根据高斯混合模型，提出了基于贝叶斯估计和高斯混合模型的小波域图像去噪算法，该方法在对自然图像去噪方法取得了较好的结果，但无法适用去SAR图像去噪中；2004年Dai等在文献“Bayesian Wavelet Shrinkage With Edge Detection for SAR ImageDespeckling.IEEE Transactions on Geoscience and Remote Sensing,2004,42(8):1642-1648”中利用加性信号相关的噪声模型，成功将高斯混合模型引入到SAR图像去噪中，并使用边缘检测算子将SAR图像进行分类，利用分类信息指导UWT域小波系数进行去噪，该方法在边缘和点目标保持方面取得了不错的效果，但对于同质区域的噪声，则无法很好的去除。

发明内容

本发明的目的在于针对SAR图像去噪中点目标难以保持和边缘容易模糊的问题，提出了一种基于变换域非局部和最小均方误差估计的SAR图像去噪方法，以在有效抑制同质区域噪声的同时，保持图像边缘和点目标清晰，提高去噪效果。

本发明的技术方案是：对含噪SAR图像进行非下采样剪切波分解，对分解得到的低频图像和高频图像分别进行PPB非局部均值滤波器去噪和MMSE估计，然后逆变换至空域得到初步去噪结果，对初步去噪结果进行分类处理，得到最终去噪结果。其实现步骤包括如下：

（1）输入一幅SAR图像Y，对SAR图像进行一层非下采样拉普拉斯分解，得到一幅低频图像Y_L和一幅高频图像Y_H；

（2）对分解得到的低频图像Y_L，利用PPB滤波器进行滤波处理，得到滤波后的图像

（3）对分解得到的高频图像Y_H进行非下采样剪切波方向滤波器分解，得到各个方向的高频子带图K为分解方向数，取值为4~18；

（4）对各个方向的高频子带图利用高斯混合模型进行建模，并利用MMSE估计进行去噪，得到去噪后的高频子带图

（5）对低频图像和高频图像进行逆剪切波变换，得到空域图像Y_Z；

（6）对空域图像Y_Z进行LCV分类处理，得到最终去噪结果Z。

本发明与现有技术相比具有如下优点：

1.本发明利用非下采样剪切波变换，可以有效的避免因变换工具缺乏平移不变性导致的去噪图像中产生的抖动失真。

2.本发明对图像低频和高频分别利用PPB滤波器和MMSE估计去噪，有效的抑制了图像在同质区域的噪声，并很好的保持了图像边缘和点目标的清晰。

仿真实验结果表明，本发明的各项客观评价指标均优于当前最好的PPB方法，在主观视觉效果方面，不但在同质区域能够很好的抑制噪声，也能保持图像边缘清晰，整体上去噪性能优越。

附图说明

图1是本发明的流程框图；

图2是本发明仿真实验使用的真实SAR图像Field；

图3是对图2标记的三个区域图像；

图4是用本发明方法与现有方法的去噪结果比较图。

具体实施方式

参照图1，对本发明的实现步骤如下：

步骤1，输入一幅SAR图像Y，对SAR图像进行一层非下采样拉普拉斯分解，得到一幅低频图像Y_L和一幅高频图像Y_H。

这里采用非下采样拉普拉斯塔型分解，是为了保持图像的冗余性质，并使分解后的图像仍具有移不变性。

步骤2，对分解得到的低频图像Y_L，利用PPB滤波器进行滤波处理，得到PPB的估计图像具体步骤如下：

2a）对低频图像Y_L的上下边界先进行M+m的行镜像扩展，再对行扩展后的图像左右边界进行M+m的列镜像扩展，得到边界扩展后的图像其中，M为搜索窗半径，m为搜索窗中所取的相似图像块的半径，M的取值范围为7~21，本实例中M取值为10，m的取值范围为1~7，本实例中m取值为3；

2b）在边界扩展后的图像中，取像素坐标范围对应在低频图像Y_L坐标范围内的任意像素点s，以该像素点为中心取一个大小为(2M+1)×(2M+1)像素的搜索窗Q^s和一个大小为(2m+1)×(2m+1)像素的图像块B^s，在搜索窗Q^s中取任意一个不同于s的像素点t为中心、大小为(2m+1)×(2m+1)像素的图像块B^t，计算图像块B^s和图像块B^t之间的相似性权值w(s,t)：

$w(s,t)=\exp [-\frac{1}{\hat{h}}\underset{r}{\mathrm{\Σ}}\log (\frac{A_r^s}{A_r^t}+\frac{A_r^t}{A_r^s})],$

其中，h为平滑因子，L为SAR图像的视数，表示以像素点s为中心的图像块B^s中的不同于像素点s的任意第r个像素点的值，表示与图像块B^s中第r个像素点坐标相对应的图像块B^t中第r个像素点的值；

2c)利用权值w(s,t)计算低频图像Y_L中像素点s的PPB估计

${\hat{Y}}_L\left(s\right)=\frac{{\mathrm{\Σ}}_{t}w(s,t){\left(A_q^s\right)}^2}{{\mathrm{\Σ}}_{t}w(s,t)},$

其中，表示以像素点s为中心的搜索窗Q^s中不同于像素点s的任意像素点q处的像素值;

2d）重复步骤2b）至步骤2c），直至遍历低频图像Y_L中所有的像素点，得到低频图像Y_L的PPB估计图像

步骤3，对步骤1得到的高频图像Y_H进行非下采样剪切波方向滤波器分解，分别得到第k方向的子带图像为k表示非下采样剪切波方向滤波器分解的不同方向，K为分解方向数，其取值范围为4~18，本发明实例中取值为18。

步骤4，对各个方向的高频子带图利用高斯混合模型进行建模，并利用MMSE估计去噪，得到去噪后的高频子带图按如下步骤进行：

4a）利用鲁棒中值法，计算第k方向子带图像的噪声标准差σ_N,k：

${\mathrm{\σ}}_{N,k}=\frac{\mathrm{median}\left(\right|Y_H^k\left|\right)}{0.6745},$

其中表示计算高频方向子代图像的中值；

4b）对第k个方向的高频子带系数利用MMSE估计去噪，得到去噪后的高频子带图

${\hat{Y}}_H^k=\frac{\left(p\left(Y_H^k\right|S=0)p(S=0)\right)({\mathrm{\σ}}_{Y_H^k}^2-{\mathrm{\σ}}_{N,k}^2)}{p\left(Y_H^k\right){\mathrm{\σ}}_{Y_H^k}^2}{Y}_H^k+\frac{\left(p\left(Y_H^k\right|S=1)p(S=1)\right)({\mathrm{\σ}}_{Y_H^k}^2-{\mathrm{\σ}}_{N,k}^2)}{p\left(Y_H^k\right){\mathrm{\σ}}_{Y_H^k}^2}{Y}_H^k,$

其中，表示高频子带图像的噪声方差，表示高频子带图像的概率密度函数，S为一个离散随机变量，p(S=0)和p(S=1)表示S取值为0和1时的概率，和为条件概率密度，其分布满足零均值高斯分布。

步骤5，对步骤（2）得到的低频图像和步骤（4）得到的K个方向高频子带图像进行逆剪切波变换，得到去噪图像Y_Z。

5a）对去噪后的K个方向高频子带做非下采样剪切波方向滤波器重构，得到一幅高频图像

5b）对步骤2所得的PPB的估计图像和步骤5a）得到的高频图像作非下采样塔型滤波器重构，得到去噪图像Y_Z。

步骤6，对去噪图像Y_Z进行分类，得到一幅最终的去噪结果图像Z，具体实现步骤如下：

6a）分别计算去噪图像Y_Z中的标准差σ和均值μ；

6b）对去噪图像Y_Z先进行F的行镜像扩展，再对行扩展后的图像左右边界进行F的列镜像扩展，得到扩展后的图像这里，F为局部图像块半径，F的理想取值范围为1~5，在本发明实例中取值为2；

6c）在边界扩展后的图像中，取像素坐标范围对应在低频图像Y_L坐标范围内的任意像素点s，以该像素点为中心取一个大小为F×F像素的局部图像块V^s，计算局部图像块V^s的局部系数方差LCV(s)：

$\mathrm{LCV}\left(s\right)=\frac{\mathrm{\σ}\left(s\right)}{\mathrm{\μ}\left(s\right)},$

其中，σ(s)和μ(s)分别为图像块V^s的标准差和均值；

6d）对去噪图像Y_Z中的像素点依据其像素值按如下准则分为两类：

当局部图像块V^s的局部系数方差LCV(s)与图像Y_Z的标准差σ和均值μ的关系满足时，定义像素点s为同质区域的点；

当局部图像块V^s的局部系数方差LCV(s)与图像Y_Z的标准差σ和均值μ的关系满足时，定义像素点s为非同质区域的点；

其中，T为收缩因子，取值范围为0~2，本发明实例中取值为0.23；

6e）重复步骤6c）~6d），直至遍历图像Y_Z中的所有像素点，则Y_Z中的所有像素点分为同质区域和非同质区域两类；

6f）在图像Y_Z中，对位于同质区域的点，取13×13像素大小的邻域进行均值滤波，对位于非同质区域的点，则保留原值不变，得到一幅最终的去噪结果图像Z。

本发明的效果可以通过以下仿真实验证明：

1、仿真条件：实验中所使用二视幅度SAR图像Field和六视强度SAR图像Bedford，这两幅图像均为256×256的真实SAR图像，图像中包含了多种场景，如农田、城镇和自然风景等，具有丰富的边缘轮廓和细节信息。

2.仿真内容：

仿真1，对上述两幅测试图，利用现有的Gamma-MAP方法对其进行去噪，去噪结果评价指标如表1和表2所示，其中，Gamma-MAP方法对图2的去噪结果如图4(a)所示；

仿真2，对上述两幅测试图，利用现有的小波域高斯混合模型滤波方法对其进行去噪，去噪结果评价指标如表1和表2所示，其中，小波域高斯混合模型滤波方法对图2的去噪结果如图4(b)所示；

仿真3，对上述两幅测试图，利用现有的非迭代PPB滤波方法对其进行去噪，去噪结果评价指标如表1和表2所示，其中，非迭代PPB滤波方法对图2的去噪结果如图4(c)所示；

仿真4，对上述两幅测试图，利用现有的迭代PPB滤波方法对其进行去噪，去噪结果评价指标如表1和表2所示，其中，迭代PPB滤波方法对图2的去噪结果如图4(d)所示；

仿真5，对上述两幅测试图，使用本发明滤波方法对其进行去噪，去噪结果评价指标如表1和表2所示，其中，本发明滤波方法对图2的去噪结果如图4(e)所示。

本发明采用的评价指标包括：

1）等效视数ENL：SAR图像中通常采用等效视数ENL来评价去噪方法的平滑效果。这里的μ₁和分别表示SAR图像同质区域的均值和方差，ENL的值越大，噪声抑制的效果越好。在本实验中，选取SAR原图中白色方形框标示的同质区域作为测试数据来计算等效视数；

2）均值：均值用于衡量去噪图像对原始图像的辐射特性的保持情况，均值越接近原始SAR图像的均值，表示辐射特性保持的越好；

3）标准差：标准差反映图像中像素的灰度值与灰度平均之间的分散程度，通过抑制性滤波，相干斑噪声的强度减弱，图像的方差减小；

（4）边缘保持指数EPI：用于考察SAR图像滤波后边缘保持效果的好坏，其定义为：

$\mathrm{EPI}=\frac{\underset{s=1}{\overset{G}{\mathrm{\Σ}}}|D_1^{\′}-D_2^{\′}|}{\underset{s=1}{\overset{G}{\mathrm{\Σ}}}|D_1-D_2|},$

其中，s表示检验像素点，G表示检验像素个数的总和，D₁和D₂分别表示原始图像中沿着边缘交界处左右或上下相邻像素的灰度值，D′₁和D′₂分别表示去噪后图像中沿着边缘交界处左右或上下相邻像素的灰度值，EPI越接近1，表示降斑方法的边缘保持性越好。在表1和表2的边缘保持指数评价的第一行和第二行，分别给出了水平方向和垂直方向边缘保持情况的指标。

表1    两幅测试图不同方法实验结果比较

表2    两幅测试图不同方法实验结果比较

从表1对SAR图像Field的实验数据可见，从本发明方法得到的去噪图像中选取图3中标记的三个区域来计算等效视数、标准差和均值，从表1中可看出，本发明三个区域的等效视数和标准差的结果均高于Gamma-MAP滤波、高斯混合模型滤波、非迭代PPB、迭代PPB滤波方法，表明本发明在同质区域噪声抑制能力要好于其他四种方法；本发明三个区域得到结果的均值比Gamma-MAP滤波稍低，但要远好于高斯混合模型滤波、非迭代PPB、迭代PPB等方法，表明本发明得到的结果和Gamma-MAP得到的结果一样，都能很好的保持原图的辐射特性；

从表2可看出，本发明方法得到的边缘保持指数的水平方向和垂直方向的指标均高于其他方法，个别与迭代PPB方法相差不大，表现出了良好的边缘保持能力。

从图2和图4可以看出，Gamma-MAP滤波器无法很好的保持图像的点目标，同质区域的噪声抑制效果也不理想，如图4(a)所显示。小波域高斯混合模型方法显示出了较好点目标保持能力，但也没能提高同质区域抑制噪声的能力，且边缘出现了较大的模糊，如图4(b)所显示。

从图4(c)可以发现非迭代PPB滤波器比小波域高斯混合模型方法显示出了更好的去噪结果，不仅能很好的保护边缘区域，而且对同质区域的平滑能力也要比小波域高斯混合模型方法效果好。但从图4(c)中还可以看出，非迭代PPB滤波器仍在一定程度上模糊了图像的边缘和点目标，对点目标的保持程度不如小波域高斯混合模型方法。

从图4(d)中可以看出，迭代PPB滤波方法要比前面各种方法去斑效果都要好，不论是在对同质区域噪声抑制能力上，还是对点目标的保护上，但该方法边缘仍会出现一定程度的模糊，且图像的均值也无法很好的保持。

从图4(e)中可以看出，本发明方法不但在等效视数方面要优于上述几种方法，且均值保持方面也体现出了较大的优势，本发明能很好的抑制同质区域的噪声，同时保持点目标的亮度不被压缩，边缘保持也很清晰。

本发明相比于现有的SAR图像去噪算法，从客观评价指标上和主观视觉质量上看都具有优越性，能够解决SAR图像去噪问题中对边缘和点目标保持困难的问题，也能很好地去除同质区域噪声，是一种可行有效的SAR图像去噪方法。

Claims (3)

1.一种基于变换域非局部和最小均方误差的SAR图像去噪方法，包括如下步骤：

(1)输入一幅SAR图像Y，对SAR图像进行一层非下采样拉普拉斯分解，得到一幅低频图像Y_L和一幅高频图像Y_H；

(2)对分解得到的低频图像Y_L，利用PPB滤波器进行滤波处理，得到滤波后的图像

2a)对低频图像Y_L上下边界先进行M+m的行镜像扩展，再对行扩展后的图像左右边界进行M+m的列镜像扩展，得到边界扩展后的图像其中，M为搜索窗半径，m为搜索窗中所取的相似图像块的半径，M的取值范围为7～21，m的取值范围为1～7；

2b)在边界扩展后的图像中，取像素坐标范围对应在低频图像Y_L坐标范围内的任意像素点s，以该像素点为中心取一个大小为(2M+1)×(2M+1)像素的搜索窗Q^s和一个大小为(2m+1)×(2m+1)像素的图像块B^s，在搜索窗Q^s中取任意一个不同于s的像素点t为中心、大小为(2m+1)×(2m+1)像素的图像块B^t，计算图像块B^s和图像块B^t之间的相似性权值w(s,t)：

$w(s,t)=\exp [-\frac{1}{\hat{h}}\underset{r}{\mathrm{\Σ}}\log (\frac{A_r^s}{A_r^t}+\frac{A_r^t}{A_r^s})],$

其中，h为平滑因子，L为SAR图像的视数，表示以像素点s为中心的图像块B^s中的不同于像素点s的任意第r个像素点的值，表示与图像块B^s中第r个像素点坐标相对应的图像块B^t中第r个像素点的值；

2c)利用权值w(s,t)计算低频图像Y_L中像素点s的PPB估计

${\hat{Y}}_L\left(s\right)=\frac{{\mathrm{\Σ}}_{t}w(s,t){\left(A_q^s\right)}^2}{{\mathrm{\Σ}}_{t}w(s,t)},$

其中，表示以像素点s为中心的搜索窗Q^s中不同于像素点s的任意像素点q处的像素值；

2d)重复步骤2b)至步骤2c)，直至遍历低频图像Y_L中所有的像素点，得到低频图像Y_L的PPB估计图像

(3)对分解得到的高频图像Y_H进行非下采样剪切波方向滤波器分解，得到各个方向的高频子带图K为分解方向数，取值为4～18；

(4)对各个方向的高频子带图利用高斯混合模型进行建模，并利用MMSE估计进行去噪，得到去噪后的高频子带图

(5)对低频图像和高频图像进行逆剪切波变换，得到空域图像Y_Z；

(6)对空域图像Y_Z进行LCV分类处理，得到最终去噪结果Z。

2.根据权利要求1所述的基于变换域非局部和最小均方误差的SAR图像去噪方法，其中步骤4所述的对各个方向的高频子带图利用高斯混合模型进行建模，并利用MMSE估计去噪，得到去噪后的高频子带图按如下步骤进行：

4a)利用鲁棒中值法，计算第k方向子带图像的噪声标准差σ_N,k：

${\mathrm{\σ}}_{N,k}=\frac{\mathrm{median}\left(\right|Y_H^k\left|\right)}{0.6745},$

4b)对第k个方向的高频子带系数利用MMSE估计去噪，得到去噪后的高频子带图

${\hat{Y}}_H^k=\frac{\left(p\left(Y_H^k\right|S=0)P(S=0)\right)({\mathrm{\σ}}_{Y_H^k}^2-{\mathrm{\σ}}_{N,k}^2)}{p\left(Y_H^k\right){\mathrm{\σ}}_{Y_H^k}^2}{Y}_H^k+\frac{\left(p\left(Y_H^k\right|S=1)p(S=1)\right)({\mathrm{\σ}}_{Y_H^k}^2-{\mathrm{\σ}}_{N,k}^2)}{p\left(Y_H^k\right){\mathrm{\σ}}_{Y_H^k}^2}{Y}_H^k,$

其中，表示高频子带图像的噪声方差，表示高频子带图像的概率密度函数，S为一个离散随机变量，p(S＝0)和p(S＝1)表示S取值为0和1时的概率，为条件概率密度，其分布满足零均值高斯分布。

3.根据权利要求1所述的基于变换域非局部和最小均方误差的SAR图像去噪方法，其中步骤6所述的对空域图像Y_Z进行LCV分类处理，按如下步骤进行：

6a)分别计算去噪图像Y_Z中的标准差σ和均值μ；

6b)对去噪图像Y_Z先进行F的行镜像扩展，再对行扩展后的图像左右边界进行F的列镜像扩展，得到扩展后的图像这里，F为局部图像块半径，F的理想取值范围为1～5；

6c)在边界扩展后的图像中，取像素坐标范围对应在低频图像Y_L坐标范围内的任意像素点s，以该像素点为中心取一个大小为F×F像素的局部图像块V^s，计算局部图像块V^s的局部系数方差LCV(s)：

$\mathrm{LCV}\left(s\right)=\frac{\mathrm{\σ}\left(s\right)}{\mathrm{\μ}\left(s\right)}$

其中，σ(s)和μ(s)分别为图像块V^s的标准差和均值；

6d)对去噪图像Y_Z中的像素点依据其像素值按如下准则分为两类：

当局部图像块V^s的局部系数方差LCV(s)与图像Y_Z的标准差σ和均值μ的关系满足时，定义像素点s为同质区域的点；

当局部图像块V^s的局部系数方差LCV(s)与图像Y_Z的标准差σ和均值μ的关系满足时，定义像素点s为非同质区域的点；

其中，T为收缩因子，取值范围为0～2；

6e)重复步骤6c)～6d)，直至遍历图像Y_Z中的所有像素点，则Y_Z中的所有像素点分为同质区域和非同质区域两类；

6f)在图像Y_Z中，对位于同质区域的点，取13×13像素大小的邻域进行均值滤波，对位于非同质区域的点，则保留原值不变，得到一幅最终的去噪结果图像Z。