CN104978716A - 一种基于线性最小均方误差估计的sar图像降噪方法 - Google Patents

Abstract

本发明公开了一种基于线性最小均方误差估计的SAR图像降噪方法。属于数字图像处理技术领域。它是一种将图像非局部相似性与稀疏表示相结合的SAR图像降噪方法。首先通过Kmeans聚类法将相似图像块聚类；再对相似块集合做奇异值分解，得到包含行列相关信息的含噪奇异值系数。为使降噪后奇异值系数更好的逼近真实系数，利用线性最小均方误差准则估计奇异值系数。接着将估计后的奇异值系数重构得到初始降噪图像块，结合初始降噪结果重新对含噪图像块进行聚类降噪，并将降噪后图像块重构得到最终降噪图像。本发明不仅去噪效果明显且能够有效的保持图像纹理细节，还具有良好的视觉效果，可用于SAR图像降噪。

Description

一种基于线性最小均方误差估计的SAR图像降噪方法

技术领域

本发明属于数字图像处理技术领域，它特别涉及基于线性最小均方误差估计的降噪方法，用于SAR图像降噪处理。

背景技术

合成孔径雷达具有全天时、全天候、可侧视成像等诸多优点而被广泛应用于多个领域，在军事打击、农林业监控、海洋开发等方面发挥着越来越大的作用。但其所成SAR图像中存在的相干斑噪声不仅影响人眼的观察还给后续的图像解译及目标识别带来了困扰，因而降噪成为SAR图像前期处理中的一个核心环节。

传统基于空域滤波方法不仅对斑点噪声的抑制能力有限而且降噪后图像的边缘细节损失严重。变换域滤波方法利用真实信号在变换域中稀疏的特性，有效的实现了变换域中图像信息与噪声的分离。这种方法的关键在于如何构造变换所对应的字典使得真实信号在变化域中更加稀疏，以及在变换域中采用那种方式将真实信息与噪声分离。最近变换域稀疏降噪与非局部相似性相结合的降噪方法取得了新突破，如BM3D方法将非局部相似图像块组合成3D图像块组，然后对3D图像块组进行3D小波变换，并使用硬阈值或维纳滤波估计真实系数，最后将系数反变换到图像域，它既利用了块内相关性又利用了块间相关性。与BM3D中相似图像块3D稀疏相似，低秩群稀疏将相似图像块组合在一起进行群稀疏表示，此时对稀疏系数进行约束不仅考虑了每个图像块本身的稀疏系数，还考虑了相似图像块间的相似系数，因此可进一步提高系数的估计精度。对于变换域中的系数，基于线性最小均方误差估计方法会获得比阈值比较方法更加精确的估计，因此被广泛应用于估计小波系数，PCA系数等。

发明内容

本发明的目在于针对现有SAR图像降噪中纹理细节损失的不足，提出一种基于线性最小均方误差估计的SAR图像降噪方法。该方法充分考虑了图像块局部稀疏与非局部稀疏的结合，将局部稀疏与非局部稀疏转换为相似图像块集合的奇异值系数稀疏，并采用线性最小均方误差估计的方法估计真实图像对应的奇异值系数，以实现在去除系数中噪声的同时更好的保护图像纹理细节所对应的小系数，最后将估计的奇异值系数重构得到降噪后的图像块集合。因此，该方法在实现SAR图像降噪的同时能更好得保持图像细节。包括以下步骤：

步骤一、相似图像块聚类

首先对含噪图像进行图像块抽取，得到图像块集合[y₁,y₂,…,y_N]，然后采用Kmeans聚类算法将相似图像块聚类，在初始降噪过程中，通过式(1)测量目标块与类中心之间的相似度：

$d(y_i,y_j)=log(\frac{z_i\left(k\right)}{z_j\left(k\right)}+\frac{z_j\left(k\right)}{z_i\left(k\right)})$       式(1)

其中y_i，y_j分别表示第i个图像块和第j个图像块，k表示图像块内的像素，在得到初始降噪图像块集合后，利用降噪后信息重新对含噪图像块进行聚类，通过式(2)测量目标块与类中心之间的相似度：

$d(y_i,y_j)=log(\frac{z_i\left(k\right)}{z_j\left(k\right)}+\frac{z_j\left(k\right)}{z_i\left(k\right)})+\mathrm{\ζ}\frac{|{\hat{x}}_i^{\′}\left(k\right)-{\hat{x}}_j^{\′}\left(k\right){|}^2}{{\hat{x}}_i^{\′}\left(k\right){\hat{x}}_j^{\′}\left(k\right)}$       式(2)

以进一步提高聚类精度。其中ζ为权重控制参数，与为初始估计的结果，d(y_i,y_j)越大表示相似度越小。

步骤二、基于局部稀疏与非局部稀疏相结合的降噪模型

图像的局部稀疏表示模型，都是基于马尔可夫随机域假设，也就是局部图像块中稀疏系数的概率密度函数与其他图像块系数相互独立。在此假设下，首先对干净图像块集合X＝[x₁,x₂,…x_m]进行局部稀疏表示，即将图像块表示成字典左乘系数的形式：

X＝D[α₁,α₂,…,α_m]       式(3)

在加性噪声模型y＝x+v下，其中v表示噪声，x，y分别表示干净图像和含噪图像，局部稀疏降噪模型等价求(4)式的最优解：

$\widehat{\mathrm{\α}}=\underset{\mathrm{\α}}{\arg \min }\left|\right|y-D\mathrm{\α}|{|}_2^2+\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}|\left|{\mathrm{\α}}_i\right|{|}_1$       式(4)

其中α_i为图像块的稀疏表示系数。为利用图像块集合中块间相似性，进一步将系数表示成字典右乘系数的形式：

[α₁,α₂,…,α_m]＝ΣΦ^T       式(5)

其中Φ^T为右乘字典，Σ＝diag{γ₁,γ₂,…γ_k}为稀疏表示系数，为对相似图像块集合同时进行局部稀疏与非局部稀疏，将图像块集合分解成式(6)的形式：

X＝DΣΦ^T        式(6)

其中D为局部稀疏字典，Φ^T为非局部稀疏字典。因此，局部稀疏与非局部稀疏相结合的降噪模型为：

$\widehat{\mathrm{\Σ}}=\underset{{\mathrm{\γ}}_i}{\mathrm{argmin}}\left|\right|Y-{\mathrm{D\Σ\Φ}}^T|{|}_2^2+{\mathrm{\λ\Σ}}_{i=1}^k{\mathrm{\γ}}_i$        式(7)

其中Y表示含噪相似图像块集合，为估计后系数。

步骤三、变换域奇异值系数线性最小均方误差估计

首先将SAR图像乘性相干斑噪声转化为加性噪声：

y＝xu＝x+x(u-1)＝x+w       式(8)

其中y，x分别表示图像的观测值和真实值，u表示相干斑噪声，由于u的均值为1，因此w为0均值加性噪声。对于相似块集合，加性噪声模型表示为：

Y＝X+W        式(9)

其中Y＝[y₁,y₂,…,y_m]，X＝[x₁,x₂,…,x_m]，W＝[w₁,w₂,…,w_m]。对Y进行SVD分解得：

Y＝DΣ_YΦ^T       式(10)

其中Σ_Y＝Σ_X+Σ_W，Σ_X与Σ_W分别为真实信号和噪声对应的奇异值，然后采用线性最小均方误差准则对真实信号的奇异值系数进行估计：

${\widehat{\mathrm{\Σ}}}_X=E\[{\mathrm{\Σ}}_X\]+Co\mathrm{\ν}({\mathrm{\Σ}}_X,{\mathrm{\Σ}}_Y)Co\mathrm{\ν}{\left({\mathrm{\Σ}}_Y\right)}^{-1}({\mathrm{\Σ}}_Y-E\[{\mathrm{\Σ}}_Y\])$      式(11)

由于相似图像块集合X和噪声W之间是不相关的，因此它们线性变换后的Σ_X，Σ_W也不相关，同时Σ_W的均值为零，式(11)可简化为：

${\widehat{\mathrm{\Σ}}}_X=E\[{\mathrm{\Σ}}_Y\]+Co\mathrm{\ν}\left({\mathrm{\Σ}}_X\right)Co\mathrm{\ν}{({\mathrm{\Σ}}_X+{\mathrm{\Σ}}_W)}^{-1}({\mathrm{\Σ}}_Y-E\[{\mathrm{\Σ}}_Y\])$         式(12)

若进一步假设Σ_X和Σ_W均为对角矩阵，则中第i个奇异值可通过式(13)进行估计：

${\widehat{\mathrm{\Σ}}}_X\left(i\right)=E\[{\mathrm{\Σ}}_Y\left(i\right)\]+\frac{{\mathrm{\σ}}_{{\mathrm{\Σ}}_X\left(i\right)}^2}{{\mathrm{\σ}}_{{\mathrm{\Σ}}_X\left(i\right)}^2+{\mathrm{\σ}}_{{\mathrm{\Σ}}_W\left(i\right)}^2}({\mathrm{\Σ}}_Y\left(i\right)-E\[{\mathrm{\Σ}}_Y\left(i\right)\])$        式(13)

其中为噪声奇异值方差。将估计系数重构得到初始降噪图像块，在初始降噪的基础上重新对含噪图像块进行聚类与降噪，并将降噪后图像块重构得到最终图像。

本发明的创新点在于在稀疏降噪过程中利用了图像的非局部冗余特性；将图像的局部与非局部同时稀疏转换为相似图像块集合奇异值系数稀疏；并利用线性最小均方误差准则估计真实图像奇异值系数，以进一步的提高估计精度，并将该方法用于SAR图像降噪。

本发明的有益效果：将图像的局部稀疏与非局部稀疏相结合，提高了稀疏降噪的性能；对相似图像块集合进行奇异值分解，得到了更适合表达相似图像块集合行和列的自适应字典，以及同时包含行列相关信息的含噪奇异值系数；并采用线性最小均方误差估计的方法估计真实图像对应的奇异值系数，能够在去除系数中噪声的同时更好的保护图像纹理细节对应的小系数，因此降噪后图像不仅去噪效果明显，同时还含有丰富的纹理细节。

本发明主要采用仿真实验的方法进行验证，所有步骤与结论都在MATLAB8.0上验证正确。

附图说明

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

图2是本发明仿真中使用的真实含噪SAR图像；

其中白色矩形区域为选择的同质区；

图3是log-KSVD方法对图2的降噪结果图；

图4是PPB方法对图2的降噪结果图；

图5是本发明方法对图2的降噪结果图。

具体实施方式

参照图1，本发明基于线性最小均方误差估计的SAR图像降噪方法，具体步骤包括如下：

步骤一、相似图像块聚类

首先对含噪图像进行图像块抽取，得到图像块集合[y₁,y₂,…,y_N]，然后从[y₁,y₂,…,y_N]中随机抽取I个图像块作为初始类中心，接着对于相似块集合中每个图像块，通过式(1)分别计算它们与I个类中心之间的距离，并将每个图像块归于类中心到自身距离最短的类，当所有图像块归类完以后，再将每类图像块集合的均值作为新的类中心，重复m次图像块归类和更新类中心过程，得到最终聚类结果。在得到初始降噪图像块集合后，将式(1)相似度测量改为式(2)，以进一步提高聚类精度。其中ζ为权重控制参数，和为初始估计的结果，同样d(y_i,y_j)越大表示相似度越小。

步骤二、基于局部稀疏与非局部稀疏相结合的降噪模型

对于每类含噪相似块集合Y，基于局部稀疏与非局部稀疏相结合的降噪模型为：

$\widehat{\mathrm{\Σ}}=\underset{{\mathrm{\γ}}_i}{\mathrm{argmin}}\left|\right|Y-{\mathrm{D\Σ\Φ}}^T|{|}_2^2+{\mathrm{\λ\Σ}}_{i=1}^k{\mathrm{\γ}}_i$      式(14)

对于式(14)模型的求解，首先对Y进行SVD分解，得到左乘字典D，右乘字典Φ^T以及含噪奇异值系数Σ，然后利用Σ估计真实奇异值系数

步骤三、变换域奇异值系数线性最小均方误差估计

首先利用式(13)估计真实图像奇异值系数，噪声奇异值方差为：

${\mathrm{\σ}}_{{\mathrm{\Σ}}_W\left(i\right)}^2=E\[{\left({\mathrm{\Σ}}_W\right(i\left)\right)}^2\]-{(E\[{\mathrm{\Σ}}_W(i)\])}^2=E\[{\left({\mathrm{\Σ}}_W\right(i\left)\right)}^2\]=\frac{{\mathrm{\σ}}_u^2}{1+{\mathrm{\σ}}_u^2}E\left(Y^2\right)$      式(15)

其中为噪声方差，L为SAR图像等效视数。根据双向方差估计理论计算真实信号奇异值方差

${\mathrm{\σ}}_{{\mathrm{\Σ}}_X\left(i\right)}^2=\max ({\mathrm{\γ}}_i^2/k-{\mathrm{\σ}}_{{\mathrm{\Σ}}_W\left(i\right)}^2,0)$       式(16)

其中γ_i表示Y的第i个奇异值，k为奇异值个数。利用式(17)将估计系数重构得到初始降噪图像块

$\hat{X}=D{\widehat{\mathrm{\Σ}}}_X{\mathrm{\Φ}}^T$       式(17)

在初始降噪的基础上重新对含噪图像块聚类并降噪，再利用式(18)将降噪后图像块重构得到最终降噪图像

$\hat{x}={\left(\underset{i=1}{\overset{N}{\Σ}}{R}_i^T{R}_i\right)}^{-1}\underset{i=1}{\overset{N}{\Σ}}\left(R_i^T{\hat{x}}_i\right)$       式(18)

其中R_i表示图像块抽取矩阵，表示降噪后图像块。

本发明的效果可以通过以下仿真实验进一步说明：

一、实验条件和内容

实验条件：实验使用的输入图像为图2，像素大小为512×512。实验中各降噪算法都使用MATLAB语言编程实现。

实验内容：在上述实验条件下，使用log-KSVD方法、PPB方法与本发明方法进行对比。降噪能力的客观评价结果用同质区等效视数ENL衡量。

实验1：用本发明方法和现有的PPB方法及log-KSVD方法分别对图2进行降噪，其中log-KSVD算法残差迭代终止条件为2.5σ²(σ²为log域转化后噪声方差)，重叠图像块大小为8×8，字典大小为64×256，训练字典迭代次数为8，降噪结果如图3所示；PPB方法的图像块大小为7×7，搜索窗大小为21×21，迭代次数为25，降噪结果如图4所示；本发明中图像块大小为8×8，权重控制参数ζ＝1,相似图像块集合数I＝60，降噪结果如图5所示。

log-KSVD算法利用稀疏表示降噪，先将含噪图像转换到log域，然后直接利用KSVD算法对其降噪，将其与PPB算法和本发明算法比较可以看出，其降噪结果不仅同质区平滑度较低，异质区纹理细节也相对模糊，整体视觉效果不佳；PPB算法利用图像的非局部相似特性降噪，同质区平滑度较高，整体降噪能力较强，但某些边缘部分仍残留噪声；本发明方法将稀疏表示与图像的非局部相似性相结合，并通过线性最小均方误差准则估计变换域稀疏系数，进一步提高其估计精度，降噪结果不仅同质区较平滑，异质区的纹理细节也得到了很好的保持，图像的整体视觉效果良好。

表1 同质区1，2的ENL指标

	含噪图像	log-KSVD方法	PPB方法	本发明方法
					同质区1	9.07	134.31	179.34	192.51
同质区2	7.91	76.08	101.83	108.30

表1给出了图2中的两个矩形同质区降噪前后的ENL情况，其中ENL值提高越多表示降噪效果越强。可以看出本发明方法略高于PPB方法，且均远高于log-KSVD方法，此结果与降噪效果图相吻合。

上述实验表明，本发明降噪方法不仅降噪效果明显，而且降噪后图像内容丰富，同时视觉效果及客观评价指标都较好，由此可见本发明对SAR图像降噪是有效的。

Claims (1)

1.一种基于线性最小均方误差估计的SAR图像降噪方法，其特征在于具体步骤如下：

步骤一、相似图像块聚类

图像中结构存在大量非局部相似信息，利用这些相似信息对图像块进行降噪不仅去噪效果明显，同时有利于保护图像的纹理细节，为了利用图像的非本地相似信息，对于含噪图像块，首先通过Kmeans聚类方法将含有相似真实信息的图像块聚类，并采用基于统计估计的相似块评价代替基于欧氏距离的相似块评价来进一步提高相似块聚类精度；

步骤二、基于局部稀疏与非局部稀疏相结合的降噪模型

为对图像块集合X＝[x₁,x₂,…,x_m]局部稀疏表示，将X表示为字典左乘系数：

X＝D[α₁,α₂,…,α_m]

其中D为字典，[α₁,α₂,…,α_m]为[x₁,x₂,…x_m]的系数集合；当X中图像块为相似图像块时，实现了X的非局部稀疏，进一步将[α₁,α₂,…,α_m]表示为字典右乘系数：

[α₁,α₂,…,α_m]＝ΣΦ^T

其中Φ^T为右乘字典，Σ＝diag{γ₁,γ₂,…γ_k}为稀疏表示系数，当同时对相似图像块集合进行局部稀疏与非局部稀疏表示时，X可表示为：

X＝DΣΦ^T

其中稀疏表示系数Σ同时包含X的行相关信息和列相关信息；因此，利用图像的非局部相似性，将局部稀疏与非局部稀疏结合的降噪模型为：

$\widehat{\mathrm{\Σ}}=\underset{{\mathrm{\γ}}_i}{\arg \min }{||Y-{\mathrm{D\Σ\Φ}}^T||}_2^2+\mathrm{\λ}{\Σ}_{i=1}^k{\mathrm{\γ}}_i$

其中Y表示含噪相似图像块集合，为估计后系数；

步骤三、变换域奇异值系数线性最小均方误差估计

为求解降噪模型，首先将SAR图像的乘性相干斑噪声转化为加性噪声，然后对含噪相似图像块集合Y进行SVD分解，对应的奇异值系数可表示为：

Σ_Y＝Σ_X+Σ_W

其中Σ_Y，Σ_X分别表示含噪奇异值系数和真实信号奇异值系数，Σ_W表示加性噪声；采用线性最小均方误差准则对真实信号的奇异值系数进行估计：

${\widehat{\mathrm{\Σ}}}_X=E\[{\mathrm{\Σ}}_X\]+Cov({\mathrm{\Σ}}_X,{\mathrm{\Σ}}_Y)Cov{\left({\mathrm{\Σ}}_Y\right)}^{-1}({\mathrm{\Σ}}_Y-E\[{\mathrm{\Σ}}_Y\])$

其中E[·]表示期望，Cov(Σ_Y)表示Σ_Y的协方差矩阵，Cov(Σ_X,Σ_Y)表示Σ_X与Σ_Y的互协方差矩阵；将估计系数重构得到初始降噪图像块，在初始降噪的基础上重新对含噪图像块聚类与降噪，并将降噪后图像块重构得到最终图像。