CN103839237A - 基于svd字典和线性最小均方误差估计的sar图像去斑方法 - Google Patents

Abstract

本发明涉及一种基于SVD字典和线性最小均方误差估计的SAR图像去斑方法：1)对输入的SAR图像中每一个图像块，根据距离公式计算该图像块与其搜寻区域内所有图像块的距离；2)根据距离最小原则构建相似集合；3)对相似集合进行奇异值分解得到SVD字典，并将相似集合向SVD字典上投影得到变换系数；4)依据线性最小均方误差准则，对变换系数进行收缩；5)对处理后的变换系数进行反变换得到去噪后的相似集合，用相似集合重构出去噪后的图像；6)迭代进行上述过程，得到最终去噪结果。本发明有效地解决了现有SAR图像去斑结果中细节保持与平滑程度的矛盾问题，在更好平滑斑点噪声同时保持SAR图像的特征信息。

Description

基于SVD字典和线性最小均方误差估计的SAR图像去斑方法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于SVD字典和线性最小均方误差估计的SAR图像去斑方法。

背景技术

合成孔径雷达（Synthetic Aperture Radar，SAR）所成的图像具有全天候、全天时、高分辨率和强大的穿透能力等特点，被广泛地应用于目标识别、变换检测和水面监视。斑点噪声是SAR成像***的一大特色，缘自基本分辨单元内地物的随机散射，在图像上表现为信号相关的小斑点，它既降低了图像的画面质量，又严重影响图像的自动分割、分类、目标检测以及其它定量专题信息的提取，因此，需要对斑点噪声进行抑制和去除。去斑方法的目标就是在去除噪声的同时保留图像的特征信息，比如保留图像的纹理、边缘和点状目标等信息。但是由于斑点噪声的乘性背景，这一目标的实现非常困难。目前，SAR图像去斑方法有多种，一般可以分为两大类，一类是基于空域的方法，一类是基于变换域的方法。

传统的空域统计类去斑方法一般先假定噪声的乘性模型，然后基于像素邻域窗口的局部统计特性进行滤波，其中具有代表性的空域滤波方法有Lee滤波、Kuan滤波、Frost滤波及其增强滤波。空域统计类去斑方法在同质区域取均值，对变化较快的点采取保留的策略，其缺点是去斑能力不足，对图像纹理的处理过于平滑或者不能有效的去除边缘周围的噪声。近年来，受非局部均值思想在自然图像去噪成功的启发，在SAR图像去斑领域涌现出了一系列基于非局部均值思想的去斑方法，其中具有代表性的去斑方法有PPB、基于贝叶斯非局部的均值滤波等。基于非局部均值思想的去斑方法，其本质上也属于空域统计类去斑方法，去斑能力非常强，但保持图像细节的的能力不够，存在易模糊图像的纹理、细节和点目标的缺点，传统的基于变换域的方法首先通过对数运算将斑点噪声的乘性模型转为加性模型，然后直接利用自然图像去噪算法去除斑点噪声。基于小波、Contourlet等变换的SAR图像去斑算法由于其多尺度、多分辨等优点获得了广泛的应用。该类方法中，由Sara Parrilli于2011年提出的SAR-BM3D方法是在PSNR指标上效果最好的方法之一。但是，变换域去斑算法本质上仍是基于固定窗口的滤波，在图像的边缘、线等区域会产生吉布斯现象。

发明内容

本发明的目的在于克服现有技术的不足，提出了一种基于SVD字典和线性最小均方误差估计的SAR图像去斑方法，以达到在增强SAR图像降斑能力的同时保持图像细节的去斑效果。

为实现上述目的，本发明所采取的技术方案:基于SVD字典和线性最小均方误差估计的SAR图像去斑方法,包括如下步骤：

1)对输入图像v进行窗口大小为3×3的Lee滤波，得到预处理后的图像v₀；

2)在图像中每隔5个像素点采样一个图像块作为中心块，根据距离公式计算各个中心块与其搜索邻域内所有图像块的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}((2L-1)\log (\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\frac{v^{\left(m\right)}\left(x_j\right)}{v^{\left(m\right)}\left(x_i\right)})+\mathrm{\γL}\frac{\left|v_0^{\left(m\right)}\left(x_i\right){-v}_0^{\left(m\right)}\left(x_j\right)\right|}{v_0^{\left(m\right)}\left(x_i\right)v_0^{\left(m\right)}\left(x_j\right)}),$

其中，L为SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素，γ是一个调和参数；

3)根据计算出的搜寻区域内所有图像块到中心块的距离，选择距离最小的20个图像块，从图像v中提取此20个图像块进行向量化形成一个含噪的相似集合Y；同样地，从预处理图像v₀中提取此20个图像块进行向量化形成一个预处理后的相似集合Y₀；

4)对Y₀中的图像块取均值，得到均值图像块y₀，对y₀取倒数得到向量β，将β作为该相似集合的权值向量；

5)利用权值向量β分别对含噪相似集合Y和预处理后的相似集合Y₀进行噪声均匀化处理，得到含有匀质噪声的相似集合

和

即

$\stackrel{\‾}{Y}={\mathrm{\β}}^{T}Y,{\stackrel{\‾}{Y}}_0={\mathrm{\β}}^T{Y}_0$

其中，T为矩阵转置；

6)将

和

分别向SVD字典上投影得到系数矩阵A和A₀，即

$A=U_0^T\stackrel{\‾}{Y},A_0=U_0^T{\stackrel{\‾}{Y}}_0$

其中，T为矩阵转置，U₀是对预处理后的相似集合

进行SVD分解得到的SVD字典，即对

进行SVD分解

${\stackrel{\‾}{Y}}_0=U_0\mathrm{\Σ}{V}^T$

得到SVD字典U₀；

7)利用系数矩阵A₀估计无噪图像在变换域的方差，即计算矩阵A₀在各行上的方差得到向量

并利用公式

${\mathrm{\σ}}_n^2\left(i\right)=\frac{D\left(n\right)}{D\left(n\right)+E^2\left(n\right)}\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\⟨d_i^2,y_k^2\⟩$

计算噪声在变换域的方差向量

其中

表示向量

的第i个元素，D(n)和E(n)分别为乘性斑点噪声的方差和均值，d_i是字典U₀的第i个原子，y_k为相似集合Y中第k个图像块对应的向量，t为相似集合的势；

8)根据线性最小均方误差准则，计算无噪相似集合在变换域的系数矩阵的估计值

即：

${\hat{A}}^i=A_0^i+\frac{{\mathrm{\σ}}_0^2\left(i\right)}{{\mathrm{\σ}}_0^2\left(i\right)+{\mathrm{\σ}}_n^2\left(i\right)}[A^i-A_0^i]$

其中，

是无噪图像在变换域的系数矩阵的估计值，表示

的第i行，

表示A₀的第i行，

表示

第i个元素；

9)利用真实系数矩阵的估计值

得到去噪后的相似集合其中，T为矩阵转置:

$\hat{Y}=y_0^T{U}_0\hat{A}$

10)按照步骤3)至9)，计算出每个含噪相似集合Y₀的无噪估计值

用各个相似集合的无噪估计值

重构出整幅图像

11)判断是否已经完成w次迭代：如果未完成，则步骤10)中的图像

取代原预处理后的SAR图像v₀，转到步骤1)；如果已经完成，则结束执行，将作为去斑结果图像。

本发明具有的优点是：

1.本发明利用的是局部自适应字典，能够很好的拟合图像的局部特性，且可以并行实现。

2.本发明由于使用了噪声均匀化操作，克服了针对匀质噪声而设计的传统的图像去噪方法运用于非匀质噪声时的天然障碍。

3.由于使用了估计真实图像块变换系数方差的方法，使得真实变换系数的估计值具有空间自适应性，同时发掘了图像的非局部自相似的特性，克服了传统方法中局部窗口估计法需要假设图像具有局部平稳性的缺点，克服了基于块的SAR图像去斑方法易过平滑图像的缺点。

4.本发明由于使用了新提出的噪声方差估计公式，能够准确地计算出SAR图像的加性非匀质噪声在字典空间的方差。

5.由于本发明对真实信号和噪声在变换域的方差估计更加准确，使得线性最小均方误差准则中的收缩系数更加准确，因而能够较好平滑斑点噪声，同时保持和恢复SAR图像的边缘和纹理细节。

以下将结合附图和实施例对本发明的技术方案进行详细说明。

附图说明

图1是本发明的流程图。

图2(a)是用于生成模拟SAR图像的自然图像Lena。

图2(b)是用于生成模拟SAR图像的自然图像Napoli。

图2(c)是单视SAR图像Rosen_SSC3。

图2(d)是六视SAR图像Toronto_MGD。

图3(a)是无噪的Lena图像。

图3(b)是PPB方法的去噪结果图。

图3(c)是SAR-BM3D方法的去噪结果图。

图3(d)是本发明方法的去噪结果图。

图4(a)是无噪的Napoli图像。

图4(b)是PPB方法的去噪结果图。

图4(c)是SAR-BM3D方法的去噪结果图。

图4(d)是本发明方法的去噪结果图。

图5(a)为单视Rosen_SSC3图像。

图5(b)为使用PPB方法对图5(a)去斑后的结果图。

图5(c)为使用SAR-BM3D方法对图5(a)去斑后的结果图。

图5(d)是使用本发明方法对图5(a)去斑后的结果图。

图6(a)为6视Toronto_MGD图像。

图6(b)为使用PPB方法对图6(a)去斑后的结果图。

图6(c)为使用SAR-BM3D方法对图6(a)去斑后的结果图。

图6(d)是使用本发明方法对图6(a)去斑后的结果图。

具体实施方式

实施例1：

本实施例参照图1，一种基于SVD字典和线性最小均方误差估计的SAR图像去斑方法，具体包括如下步骤：

1)对输入图像v进行窗口大小为3×3的Lee滤波，得到预处理后的图像v₀；

2)在图像中每隔5个像素点采样一个图像块作为中心块，根据距离公式计算各个中

心块与其搜索邻域内所有图像块的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}((2L-1)\log (\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\frac{v^{\left(m\right)}\left(x_j\right)}{v^{\left(m\right)}\left(x_i\right)})+\mathrm{\γL}\frac{\left|v_0^{\left(m\right)}\left(x_i\right){-v}_0^{\left(m\right)}\left(x_j\right)\right|}{v_0^{\left(m\right)}\left(x_i\right)v_0^{\left(m\right)}\left(x_j\right)}),$

其中，L为SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素，γ是一个调和参数；

3)根据计算出的搜寻区域内所有图像块到中心块的距离，选择距离最小的20个图像块，从图像v中提取此20个图像块进行向量化形成一个含噪的相似集合Y；同样地，从预处理图像v₀中提取此20个图像块进行向量化形成一个预处理后的相似集合Y₀；

4)对Y₀中的图像块取均值，得到均值图像块y₀，对y₀取倒数得到向量β，将β作为该相似集合的权值向量；

5)利用权值向量β分别对含噪相似集合Y和预处理后的相似集合Y₀进行噪声均匀化处理，得到含有匀质噪声的相似集合

和

即

$\stackrel{\‾}{Y}={\mathrm{\β}}^{T}Y,{\stackrel{\‾}{Y}}_0={\mathrm{\β}}^T{Y}_0$

其中，T为矩阵转置；

6)将

和

分别向SVD字典上投影得到系数矩阵A和A₀，即

$A=U_0^T\stackrel{\‾}{Y},A_0=U_0^T{\stackrel{\‾}{Y}}_0$

其中，T为矩阵转置，U₀是对预处理后的相似集合

进行奇异值分解（singular valuedecomposition，SVD）得到的SVD字典，即对

进行SVD分解

${\stackrel{\‾}{Y}}_0=U_0\mathrm{\Σ}{V}^T$

得到SVD字典U₀；

7)利用系数矩阵A₀估计无噪图像在变换域的方差，即计算矩阵A₀在各行上的方差得到向量

并利用公式

${\mathrm{\σ}}_n^2\left(i\right)=\frac{D\left(n\right)}{D\left(n\right)+E^2\left(n\right)}\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\⟨d_i^2,y_k^2\⟩$

计算噪声在变换域的方差向量

其中

表示向量的第i个元素，D(n)和E(n)分别为乘性斑点噪声的方差和均值，d_i是字典U₀的第i个原子，y_k为相似集合Y中第k个图像块对应的向量，t为相似集合的势；

8)根据线性最小均方误差准则，计算无噪相似集合在变换域的系数矩阵的估计值即：

${\hat{A}}^i=A_0^i+\frac{{\mathrm{\σ}}_0^2\left(i\right)}{{\mathrm{\σ}}_0^2\left(i\right)+{\mathrm{\σ}}_n^2\left(i\right)}[A^i-A_0^i]$

其中，

是无噪图像在变换域的系数矩阵的估计值，表示的第i行，

表示A₀的第i行，

表示第i个元素；

9)利用真实系数矩阵的估计值

得到去噪后的相似集合

$\hat{Y}=y_0^T{U}_0\hat{A}$

其中，T为矩阵转置；

10)按照步骤3)至9)，计算出每个含噪相似集合Y₀的无噪估计值

用各个相似集合的无噪估计值

重构出整幅图像

11)判断是否已经完成w次迭代，如果未完成，则步骤10)中的图像

取代原预处理后的SAR图像v₀，转到步骤1)；如果已经完成，则结束算法执行，将作为去斑结果图像。

实施例2：

本实施例中所使用的图像如图2(a)、图2(b)、图2(c)、图2(d)所示。

本实施例中，各种滤波方法都是使用matlab语言编程实现。

本实施例有两个实验组，在第一组实验中，首先向自然图像中分别加入视数L=1，2，4，16的乘性斑点噪声，然后分别使用Frost、SA-WBMMAE、MAP-S、PPB、SAR-BM3D方法和本发明实施例1的方法对它们去斑，计算去斑结果的PSNR；在第二组实验中，对真实的单视SAR图像Rosen_SSC3、六视SAR图像Toronto_MGD，分别使用PPB、SAR-BM3D方法和本发明实施例1的方法对它们去斑，比较去斑图像的视觉效果。

对于第一组实验，计算上述各种去噪方法所得图像的PSNR，如表1所示。PPB方法、SAR-BM3D方法和本发明方法的在一视模拟SAR图像上的去噪结果如图3(a)、图3(b)、图3(c)、图3(d)；图4(a)、图4(b)、图4(c)、图4(d)所示。

从表1可以明显地看出，对于各种视数的SAR图像，本发明方法和SAR-BM3D方法相比PPB方法具有绝对的优势，PSNR指标有明显的提高。比较图3(c)、图3(d)和图3(b)，明显可见PPB方法存在明显的过平滑现象，图像纹理几乎完全损失了，图像边缘等图像细节也明显变模糊，而本发明方法和SAR-BM3D方法在去除噪声的同时，能够保持图像的纹理和细节，图像的视觉效果明显优于PPB方法。仔细比较图3(c)和图3(d)可以明显地看出，SAR-BM3D方法所得图像存在块效应，图像的平滑区域明显可见板块痕迹，同时从图中人物的帽檐部分可以看出，这种板块效应也影响了图像中的边缘，而本发明方法克服了SAR-BM3D方法的这一重大缺点。

比较图4(c)、图4(d)和图4(b)，明显可见PPB方法存在明显的过平滑现象,图像纹理几乎完全损失了，图像边缘等图像细节也明显变模糊，同时丢失了大量的点目标，而本发明方法和SAR-BM3D方法在去除噪声的同时，能够保持图像的纹理和细节，图像的视觉效果明显优于PPB方法。仔细比较图4(c)和图4(d)中白框内的区域可以明显地看出，本发明方法相比SAR-BM3D具有更强的点目标保持能力。

表1各种方法去斑结果的PSNR比较

比较图5(c)、图5(d)和图5(b)，明显可见PPB方法存在明显的过平滑现象,图像纹理损失严重，同时丢失了许多点目标，而本发明方法和SAR-BM3D方法在去除噪声的同时，能够保持图像的纹理和细节，图像的视觉效果明显优于PPB方法。仔细比较图5(c)和图5(d)中的平滑区域可以明显地看出，本发明方法相比SAR-BM3D具有更强的去斑能力，同时又不会产生过平滑现象。

比较图6(c)、图6(d)和图6(b)，明显可见PPB方法存在明显的过平滑现象,图像纹理损失严重，同时在平滑区域产生了许多斑块，而本发明方法和SAR-BM3D方法在去除噪声的同时，能够保持图像的纹理和细节，图像的视觉效果明显优于PPB方法。仔细比较图6(c)和图6(d)中的平滑区域可以明显地看出，本发明方法相比SAR-BM3D具有更强的去斑能力，同时又不会产生过平滑现象。

以上实验结果表明，本发明相对于其它的一些经典的去斑方法具有更好的性能，能够更好平滑斑点噪声同时保持SAR图像的边缘、纹理细节和点目标。

上面结合附图对本发明的实施方式作了说明，但本发明并不限于上述实施方式，在本领域的普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化，均应仍归属于本发明的涵盖范围之内。

Claims (1)

1.基于SVD字典和线性最小均方误差估计的SAR图像去斑方法，其特征在于包括如下步骤：

1)对输入图像v进行窗口大小为3×3的Lee滤波，得到预处理后的图像v₀；

2)在图像中每隔5个像素点采样一个图像块作为中心块，根据距离公式计算各个中心块与其搜索邻域内所有图像块的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}((2L-1)\log (\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\frac{v^{\left(m\right)}\left(x_j\right)}{v^{\left(m\right)}\left(x_i\right)})+\mathrm{\γL}\frac{\left|v_0^{\left(m\right)}\left(x_i\right){-v}_0^{\left(m\right)}\left(x_j\right)\right|}{v_0^{\left(m\right)}\left(x_i\right)v_0^{\left(m\right)}\left(x_j\right)}),$

其中，L为SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素，γ是一个调和参数；

3)根据计算出的搜寻区域内所有图像块到中心块的距离，选择距离最小的20个图像块，从图像v中提取此20个图像块进行向量化形成一个含噪的相似集合Y；同样地，从预处理图像v₀中提取此20个图像块进行向量化形成一个预处理后的相似集合Y₀；

4)对Y₀中的图像块取均值，得到均值图像块y₀，对y₀取倒数得到向量β，将β作为该相似集合的权值向量；

5)利用权值向量β分别对含噪相似集合Y和预处理后的相似集合Y₀进行噪声均匀化处理，得到含有匀质噪声的相似集合和

即

$\stackrel{\‾}{Y}={\mathrm{\β}}^{T}Y,{\stackrel{\‾}{Y}}_0={\mathrm{\β}}^T{Y}_0$

其中，T为矩阵转置；

6)将和

分别向SVD字典上投影得到系数矩阵A和A₀，即

$A=U_0^T\stackrel{\‾}{Y},A_0=U_0^T{\stackrel{\‾}{Y}}_0$

其中，T为矩阵转置，U₀是对预处理后的相似集合

进行SVD分解得到的SVD字典，即对

进行SVD分解

${\stackrel{\‾}{Y}}_0=U_0\mathrm{\Σ}{V}^T$

得到SVD字典U₀；

7)利用系数矩阵A₀估计无噪图像在变换域的方差，即计算矩阵A₀在各行上的方差得到向量

并利用公式

${\mathrm{\σ}}_n^2\left(i\right)=\frac{D\left(n\right)}{D\left(n\right)+E^2\left(n\right)}\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\⟨d_i^2,y_k^2\⟩$

计算噪声在变换域的方差向量

其中表示向量

的第i个元素，D(n)和E(n)分别为乘性斑点噪声的方差和均值，d_i是字典U₀的第i个原子，y_k为相似集合

中第k个图像块对应的向量，t为相似集合的势；

8)根据线性最小均方误差准则，计算无噪相似集合在变换域的系数矩阵的估计值即：

${\hat{A}}^i=A_0^i+\frac{{\mathrm{\σ}}_0^2\left(i\right)}{{\mathrm{\σ}}_0^2\left(i\right)+{\mathrm{\σ}}_n^2\left(i\right)}[A^i-A_0^i]$

其中，

是无噪图像在变换域的系数矩阵的估计值，

表示

的第i行，

表示A₀的第i行，

表示

第i个元素；

9)利用真实系数矩阵的估计值得到去噪后的相似集合其中，T为矩阵转置；

$\hat{Y}=y_0^T{U}_0\hat{A}$

10)按照步骤3)至9)，计算出每个含噪相似集合Y₀的无噪估计值

用各个相似集合的无噪估计值

重构出整幅图像

11)判断是否已经完成w次迭代，如果未完成，则步骤10)中的图像

取代原预处理后的SAR图像v₀，转到步骤1)；如果已经完成，则结束执行，将

作为去斑结果图像。

Images